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1.  INTRODUCTION 


The  purpose  of  this  report  is  to  review,  unify,  and  extend  recent 
studies  of  shock  propagation  in  monatomic,  discrete  crystal  lattices. 

Most  of  these  studies,  undertaken  by  various  investigators  during  the 
past  decade  or  so,  have  employed  computer-molecular-dynamic  techniques 
to  simulate  the  motion  of  the  shock  wave  through  the  lattice,  and  this 
procedure  will  be  largely  followed  in  tlie  present  work.  Although  these 
"brute-force"  techniques  suffer  from  well-known  limitations,  the  current 
state  of  the  art  is  such  that  they  provide  the  only  practical  way  of 
treating  highly  nonlinear  problems  in  lattice  dynamics.  The  calculations 
of  the  report  will  be  confined  to  a one-dimensional  lattice.  ITiis 
limitation  has  been  introduced  partly  because  the  one-dimens ional  lattice 
is  considerably  easier  to  treat  numerically,  partly  because  we  believe 
that  most  of  the  physical  principles  are  contained  in  such  a model, 
and  partly  because  it  is  desirable  to  investigate  such  questions  as 
numerical  accuracy  before  proceeding  to  more  complicated,  three-dimen- 
sional calculations. 

The  motivation  for  studying  shock  propagation  from  a microscopic 
point  of  view  has  been  the  belief  that  a continuum  appioach  may  fail, 
in  certain  cases,  to  predict  the  effect  adequately.  Early  investigations 
of  shock  propagation  in  gases,  using  a hydrodynamic  approach,  revealed 
that  for  strong  shock  waves  the  shock  thickness  was  of  the  order  of  a 

gas  molecular  mean  free  path^.  This  condition  invalidates  the  use  of 

the  Navier-Stokes  equations^  and  prompted  investigators,  most  notably 
2 

Mott-Smith  , to  reconsider  the  shock  problem  from  the  point  of  view 
of  kinetic  theory.  More  recently,  investigations  have  been  undertaken 

using  computer-simulation  methods^. 

An  analagous  approach  for  the  case  of  shock  waves  in  solids  has 
proceeded  much  more  slowly,  owing  primarily  to  the  increased  difficulty. 
For  several  reasons,  however,  the  need  for  such  an  examination  clearly 
exists . 


1.  R.  Becker,  "Stosswelle  und  Detonation",  Z.  I’hysik  3Z1  (1922). 

2.  H.M.  Mott-Smith,  "The  Solution  of  the  Boltzmann  Equation  for  a 
Shock  Wave",  Phys.  Rev.  885  (1951). 

3.  G.A.  Bird,  "Aspects  of  the  Structure  of  Strong  Shock  Waves", 
Phys.  Fluids  13,  1172  (1970). 
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rRECKDIlG  PATS  ELANK-MOT  .-■'lU-tiD 


First,  it  is  known  that  discrete  lattices  exhibit  dispersion  which 
is  not  accounted  tor  in  the  hydrodynamic  approach.  For  example,  if  a 
one-dimensional,  harmonic  lattice  (linear  interatomic  forces)  is  subjected 
to  steady  compression,  the  resulting  "shock"  profile  does  not  approach  a 
steady  state*.  Hte  reason  is  that  the  different  normal-mode  frequencies 
present  in  the  compression  wave  travel  at  different  group  speeds;  the 
higher-frequency  modes  travel  more  slowly  than  the  lower- frequency  modes 
and  thus  trail  ever  farther  behind  the  head  of  the  front.  The  effect  is 
similar  to  the  dispersive  spreading  of  a quantum-mechanical  wave  packet. 
For  more  realistic  anharmonic  lattices  (nonlinear  interatomic  forces), 
dispersion  is  complicated  by  the  coupling  between  the  various  normal 
modes  and  it  is  important  to  determine  the  consequences  of  this  coupling. 

Second,  a f luid-dyTiamical  treatment  of  a solid  can  be  valid  only  if 
material  rigidity  can  be  neglected.  This  condition  is  probably  satisfied 
provided  the  pressure  in  the  solid  is  much  greater  than  the  shear  yield 
strength.  While  it  is  true  that  pressures  behind  strong  shock  fronts 
are  much  greater  than  the  static  yield  strength,  it  would  be  more 
appropriate  to  compare  them  with  the  yield  strength  under  dynamic 

4 

conditions.  As  Tsai  has  pointed  out,  the  yield  strength  then  is  probably 
much  greater  than  in  the  static  case,  perhaps  approaching  the  theoretical 
strength . 

Third,  continuum  treatments  of  shock  propagation  are  dependent  upon 
assumptions  about  equations  of  state  for  the  solid  under  consideration 

as  well  as  assumptions  about  the  nature  of  viscous  effects^’ The 
equations  of  state  are  inadequately  known  for  solids  and  the  origin  of 
viscous  effects  not  completely  understood.  By  use  of  computer  simula- 
tions, one  evades  the  necessity  for  these  a priori  assumptions. 


A harmonic  lattice  cannot  support  a shock  wave  in  the  usual  sense  for 
reasons  discussed  in  the  following  section.  We  shall  nevertheless 
use  this  terminology  to  refer  to  the  resulting  compression  profile 
when  the  lattice  is  subjected  to  compression. 

4.  D.H.  Tsai,  "/\n  Atomistic  Iheory  of  Sliock  Compression  of  a Perfect 

Crystalline  Solid",  in  .Accurate  Characterization  of  the  High-Pressure 
environment , edited  by  F.C.  Lloyd,  Natl.  Bur.  Stds.  Spec.  Publ. 

No.  326  (U.S.  GPO,  Washington,  DC,  1971),  p.  105. 

.5.  W.  Band,  "Studies  in  the  Theory  of  Shock  Propagation  in  Solids", 

.1.  Geophys.  Res.  695  (I960). 

6.  n.R.  Bland,  "On  Shock  Structure  in  a Solid",  .1.  Inst.  Math. 
Applications  1_,  56  (1965). 
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f-'innlly,  since  the  classic  computer  study  undertaken  by  I'ermi , Pasta, 


and  Uhuii  , considerable  speculation  has  arisen  about  the  manner  in  which 
crystal  lattices  apjn-oach  thermal  equilibrium.  The  hydrodynamic  theory 
assumes  that  thermal  equilibrium  exists  behind  the  shock  and  allows  for 
only  small  deviations  from  equilibrium  within  the  front.  It  is  clearly 
important  to  examine  the  validity  of  these  assumptions  and  they  will  be 
addressed  in  some  detail  in  this  report. 


An  attempt  to  resolve  the  points  in  the  above  discussion  appears 
to  be  of  special  importance  with  regard  to  Army-related  problems.  Tor 

g 

exiimple,  the  Zeldovich-von  Ncumann-Dor ing  (ZND)  theory  of  detonation  , 
probably  the  most  widely  used  theory  of  shock- i nduced  detonation,  is 
based  on  the  assumption  that  the  details  of  the  shock  front  can  be 
ignored.  In  the  ZNU  theory,  the  only  effect  of  the  shock  is  to  raise 
the  temperature,  density,  and  pressure  of  the  lattice  to  values  higher 
than  in  the  undisturbed  lattice.  Chemical  reactions  are  then  assumed 
to  occur  in  the  thermally  equilibrated  region  behind  the  shock  front. 
However,  if  the  shock  profile  is  not  steady,  or  if  the  approach  to 
equilibrium  is  not  much  more  rapid  than  the  times  required  for  chemical 
reaction,  the  above  assumption  must  clearly  be  called  into  question. 


Before  turning  to  the  present  calculations,  we  will  briefly  summarise 
the  existing  work  as  well  as  indicate  some  of  the  questions  that  it  has 
raised.  The  most  extensive  studies  of  shock  propagation  in  discrete  lat- 

4 h 10 

tices  have  been  undertaken  by  I'sai  and  coworkers  ’’’  over  the  last  decade. 
Thei“  more  recent  computer  codes  solve  the  atomic  equations  of  motion  for 
three-dimensional,  face-  and  body-centered  cubic  crystals  subjected  to 
shock  compression.  ;\n  essential  conclusion  of  their  work  is  that  the  shock 
profile  is  not  steady  in  time,  but  that  the  leading  edge  of  the  thermally 
equilibrated  region  behind  the  shock  propagates  more  slowly  than  does  the 
front  itself.  Consequently,  the  transition  region  grows  wider  in  time. 

The  nonstcady  behavior  has  been  interpreted  by  these  investigators  as 
being  a natural  extension  of  the  low-tcmpcrature  phenomenon  of  second 
11  12 

sound  ’ . However,  a number  of  questions  arise  regarding  this  intei’- 


7.  K.  Fermi,  J.R.  Pasta,  and  S.M.  Ulam,  "Studies  in  Nonlinear  Problems", 
l.os  Alamos  Sc i . Tab.  Rep.  TA-194fl,  1955;  also  in  Collected  Works 

of  linrico  Fermi  (IJniv.  Chicago  Press,  Chicago,  19b5),  i'.  II,  p.  97g. 

8.  B.  Lewis  and  G.  von  Elbe,  Combustion,  Flames,  and  Explosion  of  Gasc'^ 
(Academic,  New  York,  iSrSll,  Chap.  XI. 

9.  D.H.  Tsai  and  C.W.  Beckett,  "Shock  Wave  I'ropagation  in  Cubic 

I.atticcs”,  .1.  Geophys.  Res.  -bOl  (I960). 

10.  D.H.  Tsai  and  R.A.  MacDonald,  "Second  Sound  in  a Solid  Under  Shock 
Compression",  .1.  Phys.  C (^,  1,1  t1  (197.S|. 

11.  .J.C.  Ward  and  .1.  Wilks,  "Second  Sound  and  the  niermo-Mechan  i ca  1 
Effect  at  Very  Low  Temperatures",  Phil.  Mag.  4.S,  48  (1952). 

12.  M.  Chester,  "Second  Sound  in  Solids",  l’h)'s.  Rev.  1.51,  201.5  (19(v5). 


.iretatioii.  First,  current  theories  of  second  sound  suggest  that  it  i 

can  only  occur  under  extremely  limited  experimental  conditions.  1 

Specifically,  the  temperature  must  be  sufficiently  low  that  umklapp  | 

phonon  scattering  processes*'^  are  negligible,  but  sufficiently  high  | 

that  normal  processes  have  been  activated.  Unless  the  first  condition  t 

is  satisfied,  the  second-sounu  wave  is  rapidly  damped  in  time,  and,  ? 

unless  the  second  condition  rs  .met,  onl>  ballistic  propagation  occurs.  ■ 

It  is  difficult  to  see  how  such  conditions  can  be  met  in  a crystal 
under  shock  compression  where  one  would  ordinarily  expect  umklapp 

scattering  to  be  of  utmost  importance.  Second,  although  the  propagation  j 

into  the  crystal  of  a thermally  equilibrated  region  behind  the  shock 
front  is  superficially  similar  to  the  propagation  of  a second-sound 
wave,  we  view  these  processes  to  be  fundamentally  different  physical 
phenomena.  In  the  shock-wave  case,  directed  energy  which  is  deposited 
into  the  crystal  by  the  shock  front  relaxes  to  thermal  equilibrium 
behind  the  front;  if  the  relaxation  time  is  constant,  the  thermally 
equilibrated  part  of  the  crystal  will  clearly  propagate  at  the  speed 
of  the  front  and  the  profile  will  be  steady  in  time.  In  the  second- 
sound  case,  thermal  or  "random"  energy  propagates  into  the  crystal. 

Ordinarily  this  energy  propagates  only  by  diffusion  but,  under  the 
limited  experimental  conditions  discussed  above,  it  can  propagate  at 
the  higher  speed  of  second  sound.  In  our  view,  then,  in  order  to  explain 
the  nonsteady  behavior  of  the  shock  profile,  it  is  necessary  to  explain 
why  the  relaxation  time  increases  with  increasing  distance  into  the 
crystal . 

We  should  also  point  out  that  two  additional  conclusions  of 
Tsai's  work  contradict  continuum  theories  of  shock  propagation.  First, 
the  energy  density  immediately  behind  the  front  appears  to  substantially 
overshoot  its  final  equilibrated  value.  To  our  knowledge,  no  detailed 
explanation  of  this  effect  has  been  offered.  Second,  the  Hugoniot  con- 
ditions, which  relate  the  initial  and  final  values  of  the  pressure, 
density,  energy,  and  velocity  are  not  identically  satisfied.  The  lack 
of  agreement  has  been  attributed  to  the  nonsteady  behavior  and  this 
explanation  seems  reasonable  in  view  of  the  fact  that  these  conditions, 
in  their  usual  form,  are  based  on  the  assumption  of  steady  state. 

The  only  additional  computer  simulations  of  shock  propagation  in 
me  three-dimensional  lattice  with  which  we  are  aware  are  the  work  of 


l."5.  .I.M.  Ziman,  Flectrons  and  I’honons  (Oxford  University  Press, 

London,  1960),  Chap.  III. 


Paskin  and  his  coworkers  > _ These  investigations  have  been  con- 

siderably less  extensive  than  those  of  Tsai.  It  is  difficult  to 
conclude  from  their  results  whether  the  profile  is  steady  since  the 
shock  wave  has  been  allowed  to  propagate  only  a short  distance  into 
the  lattice,  i'he  Hugoniot  conditions  were  checked  and  reported  to 
be  satisfied.  Although  this  conclusion  appears  to  contradict  the  work 
of  i'sai,  we  should  point  out  that  Tsai  found  deviations  of  only  a few 
percent  and,  since  the  results  are  highly  dependent  upon  very  accurate 
determination  of  the  thermodynamic  variables  under  consideration,  the 
question  must  be  considered  unsettled. 

In  addition  to  the  work  on  three-dimensional  crystals,  shock 
propagation  in  one-dimensional  chains  has  also  been  under  active  in- 

vestigation.  Manvi,  Duvall,  and  coworkers  , in  an  effort  to 
resolve  some  of  the  more  puzzling  questions,  have  studied  the  problem 
in  detail.  They,  too,  observed  a nonsteady  profile  and  found  that  no 
temperature  rise  existed  behind  the  front.  Computer  solution  of  the 
atomic  equations  of  motion  for  a one-dimensional,  weakly  anharmonic 

19-21 

chain  have  also  been  undertaken  by  Tasi  . The  numerical  solutions 

were  then  verified  by  an  elegant  perturbation  technique.  Iasi's  result-^ 
were  in  essential  agreement  witli  tiiose  of  Manvi  et  al.  In  addition, 
the  shock  wave  was  allowed  to  propagate  considerably  farther  into  the 
lattice  than  in  previous  investigations.  At  distances  far  into  the 
lattice,  the  propagation  of  well-defined,  stable  pulses  (solitary  waves) 
was  found  in  the  vicinity  of  the  shock  front.  In  Tasi's  work  as  well 
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as  that  of  Manvi  and  his  coworkers,  it  was  assumed  that  the  atoms  of  the 
crystal  were  initially  at  rest  in  their  equilibrium  positions  prior  to 
excitation  by  the  shock  wave.  Thus,  the  seemingly  important  effect  on 
shock  propagation  of  a nonzero  ambient  temperature  was  not  examined  in 
tliese  one-dimensional  calculations. 

These  questions  will  be  addressed  in  considerably  greater  detail 
in  the  text  of  this  report  in  which  we  describe  the  results  of  calcu- 
lations undertaken  using  our  own  computer  progiam  written  in-house. 

To  summarize,  our  motivation  for  this  work  has  been  as  follows:  Shock- 

wave  studies  in  three-dimensional  crystals  have  raised,  but  in  our 
opinion  failed  to  answer,  a number  of  Impurtant  questions,  particularly 
regarding  the  existence  of  steady  state  and  the  approach  to  thermal 
equilibrium.  Attempts  to  answer  these  questions  using  simpler  one- 
dimensional  models  have  b'*en  excellent,  but  largely  inconclusive. 

In  view  of  the  importance  of  understanding  shock  propagation  to 
Army-related  problems,  it  was  considered  worthwhile  to  develop  our  own 
in-house  capability  for  studying  these  problems.  The  one-dimensional 
codi  discussed  here  is  an  initial  step  in  the  writing  of  a more  compli- 
cated th ree-d imens lona 1 code  which  we  plan  to  undertake  in  the  near 
futu re . 

ITie  organization  of  the  report  is  as  follows:  In  Sec.  2,  we  in- 

dicate the  model  as  well  as  discuss  the  interatomic  potentials  used 
in  our  calculations.  In  Sec.  ."S,  an  exact  analytic  solution  for  the 
propagation  of  compression  waves  in  the  harmonic,  one-dimensional 
lattice  is  presented.  In  Sec.  4,  the  equations  of  motion  are  rewritten 
in  nondimensional  fom  and  the  munerical  method  employed  in  the  computer 
code  discussed.  Section  5 contains  the  results  of  the  numerical 
solution  of  the  equations  of  motion  for  the  anharmonic  lattice.  Prior 
to  excitation  by  the  shock,  the  atoms  were  initially  at  rest  in  their 
equilibrium  positions.  Solitary  waves  are  found  to  propagate  in  the 
vicinity  of  the  front  and  an  analytic  approximation  for  the  shape  of 
the  solitary  waves  is  presented.  In  Sec.  fa,  the  results  of  Sec.  5 are 
extended  to  the  case  in  wliich  the  lattice  is  characterized  by  some 
nonzero  initial  temperature.  Considerable  discussion  of  the  results 
is  given  in  each  section.  Section  7 contains  a summary  and  concli.oion 
as  well  as  our  intentions  for  future  work. 


MODEL  AND  INTERATOMIC  POTENTIALS 


The  model  under  consideration  for  the  present  calculations 
consists  of  a one-dimensional  chain  of  atoms,  each  having  mass  ra, 
which  interact  through  some  interatomic  potential  to  be  specified. 
(See  Fig.  1.)  'nie  total  number  of  atoms  in  the  chain  is  N and  the 
equilibrium  spacing  between  successive  atoms  is  The  coordinate 

x.  represents  the  displacement  of  the  atom  from  its  equilibrium 

• th  . . 

atom  from  a common  origin 


J 


post  ion,  r^  is  the  distance  to  the  j 
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Figuro  1.  Model  for  simulating  shock  propagation  in  a one-dimensional, 
discrete  lattice. 


(chosen  to  be  the  equilibrium  position  of  the  first  atom),  and  r^  is 

. th  ^ 

the  corresponding  distance  to  the  equilibrium  site  of  the  j atom. 
The  coordinates  are  related  by  the  expression 


r . = r„ . + X . 
J 0]  j 


(2.1) 


The  purpose  of  the  calculation  then  is  to  solve  the  classical  equations 
of  motion  of  the  atoms  when  the  first  atom  is  subjected  to  steady 
compression  at  velocity  u. 

The  Hamiltonian  for  the  chain  can  be  written  as 


N 

H = S m V ^ <1.  (r^,  r,.  . . . rj^ 
j = l 


(2.2) 


. th 


where  vv  = dx^/dt  is  the  velocity  of  the  j particle  and  4>  is  the 
total  potential  energy  of  the  lattice.  The  differential  equation  of 
motion  of  the  particle  is  then  clearly 


d^x. 

I 

*) 

dt 


= F.  + F. 
J J 


ext 


(2.3) 


. th 


where  F.  = - Bit/Sx.  is  the  force  exerted  on  the  j particle  by  the 
^ 2 ext 

remaining  atoms  of  the  lattice  and  F.  is  the  corresponding  external 

2 ext 

force.  For  the  problem  under  consideration,  F^  is  zero  for  all 
particles  except  the  first  and  then  it  is  -F^  since  this  particle 
moves  at  constant  velocity  u. 


If  we  assume  that  the  atoms  undergo  only  small  deviations  from  their 
equilibrium  positions,  the  potential  <t  can  be  expanded  in  a Taylor 
series  about  the  atoms'  equilibrium  positions  such  that 


N 


It  = 


'’o  ^^01’  *'oj 


'2  7 ((..  . X.  x. 

^ iJ  > J 
i J = 1 


(2..1) 


where  it^^  is  a constant  which  will  arbitrarily  be  set  equal  to  zero 
hereafter  and  where 


(2.5) 


In 


The  first-derivative  term  in  the  Taylor  expansion  for  the  potential 

vanishes  since  the  expansion  is  about  the  equilibrium  positions  of  the 

atoms.  The  potential  represented  by  Eq.  (2.4)  is  the  so-called 

harmonic  potential.  Finally,  if  we  assume  only  nearest-neighbor 
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interactions  such  that  only  the  j + 1'  and  J - 1 atom  exert  an 

appreciable  force  on  the  atom,  we  have 


U 


= - Y (5; 


i . j-1 


- 2 6. 


ij 


(2.b) 


where  > is  the  force  constant  of  the  "spring"  connecting  successive 
particles  and  6 is  the  Kronecker  6.  Equations  (2.3)  - (2.b)  then  imply 

that  the  equation  of  motion  of  the  particle  is  given  in  the 
harmonic  approximation  by 


•) 

d^x  . 


dt 


^ = Y (x.  -2x,  + X.  ) + 

- ' j + 1 J j-r  J 


(2.7) 


Equation  (2.7)  is  a linear,  second-order,  differential  equation  and,  as 
we  shall  see  in  the  following  section,  it  is  amenable  to  exact  analytic 
so lut ion . 


I’articul ariy  at  low  temperatures,  the  harmonic  approximation  to 
the  interatomic  potential  is  usually  adequate  for  calculating  most 
equilibrium  properties  of  the  lattice.  For  several  reasons,  however, 
the  approximation  is  generally  inadequate  for  the  study  of  the  shock- 
wave  problem.  First,  at  the  high  temperafires  present  in  most  shock- 
wave  calculations,  the  atoms  undergo  substantial  deviations  from 
their  equilibrium  positions  and  the  iaylor  expansion  represented  in 
Eq.  (2.4)  is  no  longer  valid.  Second,  if  one  performs  a normal-mode 
analysis  to  solve  analytically  the  equations  of  motion  for  a system 
of  coupled  harmonic  oscillators,  it  is  found  that  the  energy  of  each 
mode  is  a constant  of  the  motion.  Consequently,  there  exists  no 
mechanism  in  the  harmonic  approximation  for  redistributing  the  energ>- 
among  the  various  modes.  ITiis  system  then, if  disturbed  from  thermal 
equilibrium,  can  never  return  to  equilibrium.  Finally,  it  is  known 
that  anharmonic  terms  are  responsible  for  the  steepening  effect  which 
forms  a shock  wave  from  the  initial  compression  of  the  lattice.  Any 
realistic  lattice-dynamical  study  of  shock  propagation  must,  therefore, 
be  based  upon  a model  which  includes  anharmonic  interactions. 

Two  anharmonic  potentials  which  have  been  widely  used  in  lattice- 
dynamical  calculations,  and  in  particular  in  the  study  of  shock-wave 
phenomena,  are  the  Morse  and  Lennard-.Iones  potentials,  Ihese  potentials 
are  semi-empirical  and  based  on  the  assumption  that  only  two-body 
interactions  arc  important.  If  we  further  restrict  cons i derat i on 
to  only  nearest-neighbor  interactions  the  potentials  can  be  written. 
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r 


respect ively , as 


and 


LI 


a 


6 

0 


2 


(2.8) 


(2.9) 


In  these  expressions,  c.P,  and  a are  constants  which  are  usually  fit 
to  the  experimental  data.  That  the  potentials  are  qualitatively 

similar  has  been  shown,  for  instance,  by  Choquard*"^.  Since  it  is 
slightly  easier  to  treat  numerically,  we  will  confine  our  attention  in 
the  remainder  of  this  report  to  only  the  Morse  potential. 

One  additional  potential  which  will  be  discussed  only  briefly  in 
the  text  of  this  report  is  the  so-called  hard-sphere  interaction  first 

23 

suggested  by  Northcote  and  Potts  in  their  study  of  the  Fermi-Pasta- 
Ulam  problem.  In  this  approximation,  it  is  assumed  that  the  atoms  of 
the  chain  can  be  represented  by  hard  spheres  of  diameter  d connected 
by  harmonic  springs.  Between  collisions,  the  motion  of  the  atoms  is 
governed  by  the  harmonic  approximation.  Whenever  the  atomic  separation 
of  successive  atoms  becomes  equal  to  d,  however,  the  atoms  undergo  an 
elastic  collision  in  which  their  momenta  are  interchanged . For  reasons 
discussed  in  Sec.  5,  it  is  of  interest  to  study  this  interaction  briefU. 
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3.  THi;  HARMONIC  UTTICL 


I’ropagation  of  compression  pulses  in  one-dimensional,  liarmonic 
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chains  has  been  extensively  studied  “ . For  reasons  discussed  in 
the  previous  section,  a harmonic  lattice  cannot  support  a shock  wave 
in  the  usual  sense.  Nevertheless,  it  is  interesting  to  curry  out  the 
calculation  since  the  result  lends  physical  insight  into  the  calculation 
for  the  anharmonic  case,  and  also  provides  a useful  limiting-case  check 
for  the  numerical  calculations  undertaken  subsequently. 

The  first  detailed  solution  of  this  problem  was  apparently  given 
2d 

by  Schroedinger  ' who  noted  that  the  equations  of  motion  tor  the  atoms 
could  be  put  in  a form  that  resembled  the  recursion  relation  for  Bessel 
functions  of  the  first  kind.  In  this  section  we  shall  present  a more 
detailed  calculation,  involving  a normal-mode  analysis,  which  we  believe 
is  somewhat  more  illustrative  of  the  physical  principles  involved. 

We  shall  begin  by  solving  the  atomic  equations  of  motion  for  the 
atoms  of  an  infinite,  one-dimensional,  harmonic  lattice  for  arbitrary 
initial  conditions.  Ihe  atoms  are  distributed  ssanmet ri cal ly  about  the 
zeroth  atom,  which  is  located  at  the  origin  (Sce'lig.  2).  No  external 
forces  act  on  the  chain.  We  ti'en  show  that  by  judiciously  choosing  the 
initial  velocities  and  positions  of  the  particles  for  i »»  1 , the  first 
particle  can  be  made  to  travel  to  the  right  at  a constant  velocit)’  u 
for  all  time.  In  other  woriis , the  boundary-value  problem  corres[)ond  i ng 
to  shock  proiKigation  in  a semi - i nf inite  lattice  can  bo  transformed  to 
an  i n i t i a 1 - va lue  problem  for  an  infinite  chain.  We  then  use  this  model 
to  investigate  the  shock  profile  for  two  sets  of  initial  conditions. 

In  the  first  case,  the  particles  for  which  J ^ 1 are  assumed  to  be 
initially  at  rest  in  their  equilibrium  positions,  whereas,  in  the  second 
calculation,  their  velocities  are  distributed  according  to  a Maxwellian 
distribution  .it  temperature  T. 
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convenience,  N is  chosen  to  be  odd. 


.1  General  Solution  of  the  liquations  of  Motion  for  the  Infinite, 


One-Dimensional,  Harmonic  Chain 

Consider  a chain  consisting  of  N atoms  connected  by  harmonic  springs 
of  force  constant  >.  livery  atom  has  mass  m and  is  labeled  by  index  j 
such  that 

N-1  N-1 

- — ^ J ^ -y-  • 

We  assiune  for  convenience  that  N is  odd.  It  will  be  made  arbitrarily 
large  in  the  final  results. 

rhe  differential  equation  of  motion  for  the  atom,  assuming 
only  nearest-neighbor  interactions,  is  given  by  Hq.  (2.71  without  the 
forcing  term,  viz., 

d-’x. 

m — 4 = > 2 X.  ^ x._j1  . (-S.1) 

Since  we  anticipate  oscillatory  solutions  of  this  equation,  we  assume 
a solution  of  the  form 

X.  = C.  , (.v21 

J J 

where  C.  is  the  amplitude,  uj  is  the  frequency,  and  ^ is  a phase  factor. 
The  rea-1  part  of  iq.  (.S.2)  is,  of  course,  implied.  Substitution  of 
l:q.  (.S.2)  into  (S.H  yields  for  the  eigenvectors  C.  the  relation 

- m.TC^  = > (C.^j  - 2 C.  . C._j)  . (.S..S1 

Considerable  simplification  in  the  solution  of  I'.q.  (..S.S1  arises  if 
we  .issume  periodic  boundary  conditions  such  that 

C.  = t:.  . (->.4) 

.1  .1 

As  N is  made  arbitrarily  large,  the  boundary  conditions  at  the  ends  of 
the  chain  are  insignificant.  I'o  solve  liq.  (.S..^1,  let 


C.  - zJ  , 


(S.S) 


which  upon  substitution  into  f.q.  (.^.11  yields 


1 


27Tis'/N  , -N+1  -N+3  N-1 

I = e , s'  = — — • — — • • • • ~T 


(3.6) 


isa,, 
z = e 0 


Here , 


(3.8) 


where  L=Na^  is  the  total  length  of  the  chain  and  is  the  equilibrium 
spacing.  Finally,  from  Fqs.  (3.7)  and  (3.5)  we  have 


d = e'J"% 


The  superscript  s has  been  added  to  denote  that  Cj  is  the  eigenvector 
associated  with  a particular  one  of  the  N values  of  s. 

If  the  solution  represented  by  Fiq.  (3.9)  is  substituted  into  the 
original  expression,  F.q.  (3.3),  we  obtain  after  some  algebra  the  result 


where 


= s 4>7m^ . 


(3. 10) 


(3.11) 


Ilie  frequencies  represented  by  Fq.  (3.10)  are  the  N normal-mode 
frequencies  in  which  the  lattice  can  oscillate.  The  speed  v^  with 

which  a particular  normal  mode,  characterized  by  frequency  , propa- 
gates is  given  by  the  relation 


3cOs  ‘"O^'O 


|cos(saQ/2) 


(3.12) 


The  low-frequency,  long-wavelength  modes,  characterized  by  small  values 
of  s,  propagate  at  a velocity  nearly  equal  to  (jj^^aj^/2  which  is  just  the 

ordinary  sound  speed  in  the  crystal.  Hie  higher-frequency  modes,  however, 
propagate  more  slowly  and,  thus,  the  lattice  exhibits  dispersion. 


The  solution  represented  by  t'q.  (3.2)  is  only  a particular 
solution  for  x..  In  general,  in  order  to  satisfy  the  initial  con- 
ditions, must  be  expanded  in  terms  of  the  complete  set  of  eigen- 
vectors . Thus  we  have 
J 


s * *s  (3.13) 

. e + n C . e 

J s j J 


where  the  real  part  has  been  taken  explicitly  and  where  the  asterisks 
denote  complex  conjugates.  The  quantities  n^  and  6^  are  to  be  deter- 
mined from  the  initial  conditions. 

l.et  CL.  and  a,  represent  the  values  of  x.  and  dx./dt  at  time  t=0. 

J J J J 

Substitution  into  T!q.  (3.13)  then  yields 

" "s  S*' 

and  ' (3.14) 

= 2 2^  ['’s  “s  ‘ j ' ^s  "'s  S J • 

*s 

If  Hqs.  (3.14)  are  multiplied  by  C.‘  and  summed  over  j,  one  finds  with 
the  help  of  the  orthogonality  relation 

Vt;^  = N6  , (3.15) 

^ J J ss' 


(3.15) 


n,  4 E f j - ■ 


Substitution  of  Hq.  (3.1b)  into  Hq.  (3.1.3)  then  yields 


(3. lb) 


^ ^ { '^k  ^ ‘^s 

s , k 

) 

^ _ sin  [(j-k).saj,  + ui^t  ] J . 


1.3.  r) 


p 


In  obtaining  this  result,  we  have  used  Eq.  (3.9). 

Equation  (3.17)  represents  the  general  solution  of  the  equation 

of  motion  for  the  displacement  of  the  atom  in  terms  of  the  initial 
conditions  of  all  the  particles  in  the  chain.  The  velocity  can,  of 
course,  be  obtained  by  differentiation  with  respect  to  t.  As  N becomes 
arbitrarily  large,  the  right-hand  side  of  Eq.  (3.17)  changes  negligibly 
as  s assumes  successive  values  and  the  sum  over  s can  be  replaced  by 
an  integral.  Examination  of  Eqs.  (3.6)  and  (3.8)  shows  that  this 
replacement  is  effected  by  using  the  prescription 


ITie  region  of  s space  contained  between  ± is  known  as  the  first 

Brillouin  zone.  I'hus  Eq.  (3.17)  becomes 


(3.18) 


It  remains  only  to  evaluate  the  integral  appearing  in  Eq.  (3.18). 
To  do  so,  we  make  the  change  of  variable 


sa 


y 


0 


Some  rearrangement  of  E.q.  (3.18)  then  yields  the  equivalent  expression 
n/. 


X:  =e  E 


J IT 


J dy  cos  [2(j-k)y]  cos(w^^t  sin  y) 


(3. 19) 


cos  [2(j-k)y]  sin(u,Qt  sin  y)  | 


2J 


in  terms  of  y.  We  now  make  use 


We  have  used  hq.  (.5.10)  to  express  uj 
ot  the  generating  functions 


cos  (z  sin  O')  = 0 [z)  + 2 ^ cos  (2k6) 

^ k=l 


sin  (z  sine)  = 2 J,j^^^(z)sin  [(2k+l)eJ 


(3.20) 


where  .)  is  the  Bessel  function  of  the  first  kind  of  order  H.  Use  of 

i 

Kq.  (3.20)  in  (3.19)  allows  immediate  evaluation  of  the  y integration. 
Ultimately,  the  result  becomes 

^.'0 

= Zi  K •’2j-2k^"o'^ 

(3.21) 

>V>]  ■ 

m= j -k 

Tiie  velocity  of  the  particle  can  be  obtained  by  time  differentiat  ion 
of  l;q.  (3.21).  Noting  that'^*^ 


"n>=)  r 1 

-IE-  ■ [''..-1'  = ' - 

allows  us,  after  some  manipulation,  to  write 

oo 

K=-od 


(3.22) 

the  result  in  the  form 


J3.23) 


■^2j-2k-l^‘^0^^ 


liquations  (3.21)and  (3.23),  which  represent  a general  solution  to 

the  equations  of  motion  for  the  particle  in  the  harmonic  chain,  .ire 
the  central  result  of  this  subsection.  In  the  following  subsection 
they  are  applied  to  the  shock-wave  problem. 
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Application  to  Shock-Wave  Problem 

In  order  to  apply  the  solution  for  the  infinite  chain  to  the 
shock-wave  case,  we  must  generate  the  following  set  of  initial  condi- 
tions: (1)  The  initial  displacements  and  velocities  to  the  left 

of  the  first  particle,  which  are  at  our  disposal,  must  be  chosen  so  that 
the  first  particle  moves  at  constant  speed  u for  all  t.  C2)  The 
initial  displacement  cf  the  first  particle  is  zero  and  its  initial 
velocity  is  u.  (3)  The  initial  conditions  for  particles  to  the 
right  of  the  first  particle  are  arbitrary.  Under  these  conditions, 
the  steady  compression  of  the  first  particle  will  generate  a shock 
wave  in  the  right-most  part  of  the  lattice. 

The  symmetry  of  the  problem  suggests  that  we  examine  the  follow- 
ing initial  conditions  for  the  particles  to  the  left  of  the  first 
particle  in  the  chain 


a.  = - a , T 
k -k+2 


a,  = 2u  - a , T 
k -k+2 


(3.24) 


Note  that  these  conditions  imply  = 0,  = u.  If  these  initial 

conditions  are  substituted  into  Eqs.  (3.21)  and  (3.23),  we  find  after 
some  tedious  manipulation  of  the  summation  indices 


CO 

k=0 

k=2  m=j -k 

oo 

^ ^k['^2j-2k  " '^2j+2k-4  ] 


and 


i 


3 

i 

i 


V . ( t ) = u .1  , . , + 

J Jj-2 


k=: 


[■^2i-2k^‘^()^'  ■ ''2j^2k-4^“0^^] 


C3.2b) 


0 ^ r I f 

■ 2 ^k  L‘’2j-2k-<-l '■‘^0"'  ''2j-2k-r'“0 


,t)  - .1, 


•’2j  + 2k-5*^‘^u’^^  "■  '^2jH-2k-5'-“0^^ 


For  the  first  particle  in  the  chain,  v;e  substitute  j=l  in  Fqs.  (o.25) 
ami  (o.2b)  anJ  obtain  the  result 


Xj I t ) = ut 


Vj(tl  = u 


(0.27) 

(0.28) 


for  all  t.  In  obtaining  liqs . (.^.27)  and  (.'5.28)  we  have  made  use  of  the 


.^0 


following  relations  for  Bessel  functions'  : 


.1  (X)  = (-1)  .1  (X) 

-n  ' n 


.J^(x, 


k = l 


(.1.29) 

(.8..i0) 


and 


E 


k=() 


(2k+l)  d,^^,(x)  = X 


(.1.  .81  ! 


Thus,  the  initial  conditions  suggested  in  Fq.  (.8.24)  give  rise  to 
steady  motion  of  the  first  particle  to  the  right  at  velocity  u. 
liquations  (.8.2.S)  and  (.8.2b)  represent  the  time-dependent  velocity  and 


. th 


displacement  of  the  j particle  in  response  to  the  resulting  shock 
wave.  Our  interest,  of  course,  is  only  in  particles  to  the  right  of 
j=0. 


3 . 3 Propagation  in  the  Initially  Quiescent  Lattice 


Kquations  (3.25)  and  (3.26)  can  be  greatly  simplified  if  we  assume 
that  initially  all  particles  to  the  right  of  the  first  are  at  rest  in 
their  equilibrium  positions.  Linder  this  assumption,  tne  equations 

become  simply 


2u 


(2k+l)  J2j  + 2k-l^“'0*^^ 

k=0 


and 


Vj(t)  u J2j_2^“'0 


mpt)  + 2u  X]  ‘J2j+2k-4^V^ 
k=2 


(3.32) 


(3.33) 


The  infinite  series  in  Kqs . (3.32)  and  (3.33)  converge  rapidly  and 
evaluation  of  the  sums  poses  no  major  difficulty.  We  have  written 
a short  computer  program  which  performs  the  calculation  and  we  now 
turn  to  a discussion  of  the  results. 


In  Fig.  3,  we  have  plotted  the  nondimensionalized  velocity  of  each 
particle  given  by 

Vj  = Vj/u  (3.34) 

as  a function  of  particle  number,  j.  Results  are  presented  for  two 
different  values  of  the  nondimensionalized  time 


T = Uyt  . (3.35) 

From  the  results  it  is  clear  that  the  shock  front  propagates  at  the 
rate  of  approximately  one-half  particle  per  unit  of  nondimensional 
time  T.  This  result  could  have  been  anticipated  from  Eq.  (3.12),  which 
indicates  that  the  maximum  normal-mode  velocity  is  one-half  lattice 
spacing  per  unit  of  nondimensional  time.  The  most  important  conclusion 
that  can  be  drawn  from  the  figure  is  that  the  shock  profile  is  not 
steady  in  time.  Rather,  behind  the  front,  there  exists  a region  of 
oscillation  with  each  particle  oscillating  about  a mean  value  of  u,  the 
compression  velocity.  The  extent  of  the  region  of  oscillation  increases 
with  increasing  distance  into  the  lattice  and  is  indicative  of  the 
dispersive  character  of  the  lattice  discussed  previously.  Since  the 
higher-frequency  components  of  the  compression  will  continue  to  trail 
farther  behind  the  head  of  the  front,  we  conclude  that  the  region  of 
oscillations  must  continue  to  grow. 

In  Fig.  4,  we  have  shown  the  velocity-time  trajectories  of  two 
typical  particles  in  the  lattice.  Prior  to  excitation  by  the  shock 

front,  the  100^*'  particle  is  at  r»'St  in  its  eijui  librium  position.  Upon 


imerisional  velocity  as 


I 


excitation,  it  reaches  a maximum  nondimens ional  velocity  of  about  1.25, 
and  then  oscillates  with  decaying  amplitude  about  unity.  A similar 

effect  is  observed  for  the  SOo'"^  particle  in  the  lattice,  except  that 

the  decay  time  for  the  velocity  is  longer  than  for  the  100^*^  particle, 
again  because  of  dispersion. 


Well  behind  the  shock  front,  all  particles  are  observed  to  move 
uniformly  at  the  compression  velocity;  consequently  there  is  no  tem- 
perature rise  behind  the  shock  front.  This  effect  is  not  unejcpected 
because,  as  pointed  out  earlier,  there  is  no  mechanism  in  the  harmonic 
lattice  for  "randomizing"  the  energy  deposited  in  the  lattice  by  the 
shock  wave.  That  the  velocities  behind  the  front  should  eventually 
approach  the  compression  velocity  can  be  seen  directly  from  the 
asvTnptotic  expansion  (large  t)  of  Fq.  (3.33).  For  v fixed  and  finite, 
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one  has 

tim  J (z)  = 0 . (3.36) 

V 

2-HX> 

In  Fq.  (3.33),  let  j+k=?..  Some  rearrangement  of  the  equation  then 
yields  the  equivalent  expression 


v.(t)  = u .I.j.,(u.yt)  . 2u  2 X] 

^ 1=7, 


(3.37) 


Fquations  (3.30)  and  (3.36)  then  imply 

Him  Vj (t)  = u (3. 38) 

•t-x» 

j finite 

3 . 4 Fffect  of  Nonzero  Initial  Conditions 

In  the  preceding  subsection  we  have  obtained  a solution  of  the 
atomic  equations  of  motion  for  the  shock-compressed  chain  under  the 
assumption  that  the  atoms  were  initially  at  rest  in  their  equilibrium 
positions.  Different  microscopic  initial  conditions  will  lead  to 
different  results,  of  course,  and  it  is  of  interest  to  ask  if  an 
average  taken  over  a large  set  of  such  initial  conditions  will  lead 
to  a steady  profile. 

To  examine  this  question,  consider  an  ensemble  of  initial  conditions 
in  which  the  particles  in  the  right-most  part  of  the  chain  are  in  their 
equilibrium  positions,  but  their  velocities  arc  characterized  by  a 
Maxwellian  distribution  at  temperature  T.  These  initial  conditions  are 
summarized  by  the  initial  distribution  function 
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F(0)  = - n exp  - m a . / (2kT)  J 
Z j=2 


(.^.39) 


where  k is  Boltzmann's  constant  and  Z is  the  partition  function  which 
normalizes  the  distribution  to  unity.  The  ensemble  average  of  x.(t) 

and  v.(t)  is  then  obtained  by  multiplying  liqs.  (3.25)  and  (3.26)  (with 

all  set  equal  to  zero)  by  Kq.  (3.39)  and  evaluating  the  integral  over 

phase  space. 

For  the  velocity,  we  obtain  the  result 


< v,(t)  > 


ao  <x> 

"i  /“-’  / ‘'S  ' 

_ JD  ^ 


2k--l^“0^’  * ^k  '’2j-2k^'^n'^ 


(3.40) 


‘^k  •^2j*-2k-4^“o' 


Ci 


ex})  [-m  (X  / ( 2ki  ) ] 


where  the  brackets  denote  ensemble  average.  The  integrals  over  the 
vanish  because  of  the  antisymmetry  of  the  integrand,  leaving 


< V . (t  1 
J 


> = Tu  d_,._,(^j,t)  > 2u  ^ .1,.^ 

i = •> 


2k-4^V^ 


(XI 

oo  ^ 

yj  da,  f d a....  exp  [-  m a.  /(2kT)]  . 


(3.41 ) 


The  remaining  integration  in  Fq.  (3.41)  is  just  the  definition  ol  the 
partition  function.  Consequently,  we  again  obtain  the  result  of 
F(|.  (3.33).  A similar  result  holds  for  the  ensemble  average  of  the 
displacement,  Xj(t). 


We  conclude  that  an  average  taken  over  many  experiments  in  which 
the  initial  conditions  are  characterized  by  hq.  (3.39)  will  lead  to 
the  same  results  for  the  (average)  velocity  and  displacement  as  for 
the  case  in  which  the  initial  conditions  are  zero.  In  essence,  averaging 
the  initial  conditions  and  then  doing  an  experiment  is  equivalent  to 
averaging  the  results  of  many  experiments  with  different  initial 
conditions.  This  is  certainly  not  true  in  general  and  results  from  the 
fact  that  the  solutions  of  the  ecjuations  of  motion  depend  only  linearl)' 
on  the  initial  values  of  the  velocity  and  displacement. 

4.  NONDIMKNSlONALIZIiD  EQUATIONS  AND 
METHOD  FOR  SOLUTION 

The  equations  of  motion  for  the  atoms  of  a one-dimensional, 
harmonic  lattice  can  be  solved  analytically,  as  has  been  shown  in  the 
last  section.  For  the  case  in  which  anharmonic  terms  are  present,  how- 
ever, one  must  resort  to  numerical  techniques.  In  thi^l^sect ion , we 
consider  the  differential  equation  of  motion  for  the  j particle  when 
the  interatomic  potential  is  the  Morse  potential  represented  by  F,q.  (2.8). 
We  define  certain  nondimens ional  quantities  and  write  the  equations  in 
terms  of  these  quantities.  The  numerical  technique  used  in  the  computi-r 
solution  of  these  equations  is  then  discussed.  Finally,  the  mechanics 
of  the  computer  code  and  the  actual  quantities  calculated  are  indicated 
briefly. 

Consider  the  Morse  potential  in  F.cj.  (2.8)  which  represents  the 
interatomic  interaction  between  nearest-neightbor  jjarticles  in  the 
lattice.  As  the  values  of  x.  become  small,  the  exponential  term  can 
be  expanded  to  produce,  to  lowest  order, 

’ A 

4 --  [)  V (X.  - X )“  , (4.11 

m A"  1 1 - 1 ^ 

1 = j 

which  is,  of  course,  just  the  harmonic  limit  of  this  i nter.iction . We 

■) 

can,  therefore,  equate  the  constant  |)a“  to  one-half  the  harmonic  tT)rce 
constant,  y.  Use  of  F!q . (3.11)  then  implies 


where  is,  again,  the  maximum  normal -mode  frequency  of  the  correspondi ng 


harmonic  lattice. 

The  total  potential  energy  of  the 

l.'ittice  c.in  then 

be  written 

.S.N 

■) 

"'“o'’  'ST'  1 » > 

= -T  2^  [- 

“ 8a  ’ 


where  we  have  used  1;l(S.  (d.8)  and  (4.2).  In  the  abs^jice  of  external 
forces,  the  differential  etjuation  of  motion  of  the  j particle  is  just 
given  by  the  result 

2 2 
d"x.  -34  mu)Q 

=4^  {‘-‘"P  [*  ‘•'"P 

dt  J ' 

- exp  [-  2a(x.^j-x.)]  + exp  [-a  (x  . ^ ^ -x  . )]  | (4.4) 

where  the  right-hand  side  has  been  obtained  by  differentiation  of  hij. 
(4.5)  with  respect  to  x^ . Since  the  above  equation  contains  no  external 

forces,  it  applies  only  for  j **  2.  The  first  particle  which  moves  at 
constant  velocity  u for  all  time  satisfies  the  equation  of  motion 


with  initial  conditions  Xj  = 0,  ^^j/dt  = u. 

liquations  (4.4)  and  (4.5)  constitute  a set  of  N, coupled,  nonlinear, 
second-order  differential  equations  which  must  be  solved  numerically 
for  various  values  of  the  parameters  under  consideration.  Before 
discussing  the  method  of  solution  of  these  equations,  it  is  most  con- 
venient to  convert  them  to  a set  of  first-order  equations  and  to  rewrite 
them  in  terms  of  certain  nondimens  iona^l^qiiant  i t i es . To  do  so,  we  define 
a nond imens iona 1 displacement  of  the  j particle  by 


.s  . = — X . , 

J i‘  .1 


and  a nond imens iona 1 time,  r,  by  the  relation 

t = w^jt.  (4 

The  nond  imens  ional  velocity  of  the  jiarticle  is  then  given  by 

d.S^  J dx . 

\i  dt  u dt 

or  simply  the  real  velocity  normalized  by  the  compression  velocity. 
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It  wl'  make  use  of  fqs.  (4.b) 
in  the  form 


(4.8),  we  can  cast  liqs. 


(4.4)  and  (4.5) 


S.  = V. 
J 


V . 


.1 


IT-  I [-  - s )]  - e.xp  [-  A_^(ii  - S )] 

ml  ■ - 

- e,xp  [-  - S.)]  - exp  [-A^(S.^j  ' j ^ ‘ 


V 


1 


0 


(4.9) 


In  these  ecjuations,  is  given  by  the  expression 


A 

m 


gu 

10 

o 


14.10) 


and  each  dot  represents  differentiation  with  respect  to  the  dimension- 
less time  T.  The  initial  conditions  of  the  first  particle  are  given  h\' 

Sj  = 0,  V'j  = 1 and,  for  the  remaining  particles,  the  initial  conditions 

depend  upon  the  initial  state  of  the  lattice  prior  to  compression.  Note 
that  the  original  N second-order  equations  have  been  converted  to  JN 
first-order  eipiations. 

Most  investigators  have  not  employed  this  method  of  nondimens ion- 
alizing  the  equations  of  motion,  but  have  used  the  lattice  constant, 
to  normalize  the  displacement.  As  a result,  it  is  necessary  that 

they  supply  two  parameters  as  input  data  for  their  numerical  calculations. 
The  method  which  we  have  used  requires  spec i f ica t i on  of  onl\-  one  parameter, 
namely,  A^^,  to  solve  the  equations  and  appears  more  convenient  to  use 

in  the  numerical  ca 1 cu 1 at i ons . It  is  interesting  to  point  out  that 

Manvi  et  al.*  found  it  rather  surprising  that  the  solution  of  the 
equations  of  motion  appeared  to  depend  only  oti  the  product  of  the 
parameters  g and  u.  The  normalization  adopted  above,  bowever,  reveal', 
that  this  must  clearly  be  the  case. 

To  solve  bq.  (4.9)  we  employed  a fourth-order  Runge-kutta  scheme'^' . 
iliven  the  values  of  the  functions  on  the  left-hand  side  of  l.q.  (4.9) 
at  time  t,  this  method  a]iproximates  their  values  at  time  t ♦ .'.t  by  .i 

.51.  . R.  Carnahan,  H.A.  l.uther,  and  I.O.  Wilkes,  i ed  Numerical 

Methods  (Wiley,  New  York,  19('9),  Chap.  ii. 


fourth-order  polynomial  in  At.  The  2N  equations  in  Eq.  (4.9)  can  be 
written  in  the  general  form 


"ri  = ^2’  •••  ^ • (4.11) 


If  y.  . represents  the  value  of  the  function  y.  at  the  beginning  of  the 


th 


' J.t 


J 


time  interval,  its  value  at  the  end  of  the  interval,  y 


by  the  following  algorithm 
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j.i^l’ 


IS  given 


where 


■'^j.i+l 

= S.i 

+ At(k.,  + 2 k.^  + 2 k.,  + k.  )/6 
jl  j2  j3  j4' 

(4.12) 

Si  = 

^'2,i"  • ■^'2N,i^ 

(4.12a) 

1/2  Ar  k.^ 

(4.12b) 

N:  ■ 

’ ^ 2,i" • 2N,i^ 

(4.12c) 

V . . = 

' J . 1 

V.  . -»• 

' J ,1 

1/2  At  k.. 

(4.12d) 

k.,  = 

^2,i‘  • "Sn,!^ 

(4.12e) 

_ * 

' j .1 

= V . + 

' J . t 

At  k . , 

(4. 12f) 

S-1  = 

’ 2,i"  • 2N,i^ 

(4.12g) 

i'hc  initial  values  of  all  functions  are > of  course,  known  and  successive 
application  of  Eqs . (4.12)  allows  us  to  calculate  them  at  any  time  t. 


We  should  mention  that  in  our  initial  calculations,  we  did  not  use 
the  Kunge-Kutta  method  discussed  above,  but  rather  the  improved  Euler- 


Cauchy  method" 
the  value  y.  . , 

j , 1 + 1 

the  value  of  the 
evaluated  ;ind  an 


used  by  Tsai  et  al.  and  by  Manvi  et  al.  In  this  method, 
is  approximated  by  a linear  function  of  At,  based  on 


function  and  its  slope  at 
improved  value  of  v.  . , 

' j , 1 + 1 


i.  The  slope  at  i+1  is  then 
calculated  from  the  average 


of  the  slopes  at  i and  i+1.  The  process  is  then  repeated  until  success- 
iv'e  iterations  agree  to  within  some  prescribed  tolerance.  It  is  known, 
however,  that  this  iterative  technique  is  not  so  accurate  as  the 
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fourth-order  Runge-Kutta  method  for  the  same  value  of  At.  In  practice, 
we  found  that  although  a number  of  checks  used  to  validate  the  computer 
code  (discussed  later  in  this  section)  were  satisfied,  varying  the  step 
size.  At,  led  to  completely  different  profiles  when  the  Euler-Cauchy 
method  was  used.  The  effect  became  more  and  more  pronounced  as  both  t 
and  the  nonlinearity  parameter,  A , increased.  Some  further  discussion 
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of  this  fact  has  been  given  by  M.mvi 


In  addition  to  the  velocities  and  displacements  of  each  of  the 
particles  in  the  chain  as  a function  of  time,  we  have  calculated  in 
the  program,  a number  of  other  quantities  of  interest  such  as  the 
kinetic  and  potential  energies  of  each  particle  and  the  force  exerted 


on  the  particle  by  the  N+l'^*’, 


All  energies  have  been  nondimen- 


sionalized  by  the  kinetic  energy  of  the  first  particle,  1/2  mu“  , and 
all  forces  by  the  factor,  1/4  muiyU.  The  potential  energy  of  a single 


particle  was  defined  simply  to  be  one-half  of  the  potential  energy  of 
each  of  the  "springs"  to  which  the  particle  is  attached.  Dius , from 

Eq.  (4.3),  we  have  for  the  potential  energy  of  the  i*^^  particle. 


(‘1.  13) 


This  is  not  the  only  possible  operational  definition,  but  is  a reason- 
able one. 

In  the  remainder  of  this  report,  unless  otherwise  indicated,  when  wo 
refer  to  the  quantities  discussed  above  we  shall  mean  their  nondimen- 
sionalized  values.  For  reference.  Table  I indicates  the  quantity  under 
consideration  (column  1),  the  s>Tnbol  for  its  nondimensional  value 
(column  2),  the  normalizing  factor  by  which  the  real  quantity  is 
divided  to  make  it  nondimensional  (column  3),  and  the  appropriate 
formula  for  the  quantity  of  interest,  where  applicable  (column  4). 


32.  A.  Ralston,  A First  Course  in  Numerical  Analysis  (Metiraw-Hi 1 1 , 

New  York,  U)b5),  Chap.  S. 

33.  R.  Manvi,  "Shock  Wave  Propagation  in  a Dissipating  Lattice  Mc'del", 
I’h.l).  Thesis  (Washington  State  University,  1908)  (unpublished). 


TABLE  I.  NO.NDIMENSIONAL  PARAMETERS  IN  COMPUTER  CALCULATION 


The  remainder  of  this  report  will  be  concerned  primarily  with  the 
numerical  solution  of  Eq.  (4.9)  using  the  method  just  outlined.  The 
results  are  discussed  in  the  following  two  sections.  We  conclude  this 
section  by  indicating  a number  of  checks  we  used  to  insure  the  accuracy 
of  our  computational  results. 


First,  we  took  the  harmonic  limit  of  Eq.  (4.9)  by  expanding  the 
exponentials  and  retaining  only  the  first  nonvanishing  terms.  The 
equations  then  become 
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(4.14) 


with  the  Siime  initial  conditions  as  before.  This  set  of  equations  was 
solved  for  the  initially  quiescent  lattice,  using  a step  size  of 
At  = 0.05,  and  produced  the  same  numerical  results  as  obtained  for  the 
analytic  solution  of  Sec.  3.3.  A similar  result  can  also  be  obtained 
by  simply  making  very  small  (say,  0.001)  in  Eq.  (4.9). 


Sicond,  in  all  our  calculations  we  calculated  the  work  done  by  the 
external  force  needed  to  move  the  first  particle  at  a constant  velocity 
of  unity.  This  external  force  is,  of  course,  just  the  negative  of  the 
force  exerted  on  the  first  particle  in  the  lattice  by  the  second.  A 
simple  calculation  reveals  that  in  an  infinitesimal  time  interval  of 
At,  the  external  force  does  an  amount  of  work  W,  normalized  by 
*) 

1/2  mu'",  given  by 


. exp  [-  2A^(S,^  - Sj^j  - exp  A^(S,^.  - S^^)] 

X (S„-Sj.) 

where  i and  f denote  the  values  of  Sj  and  S.,  at  the  beginning  and  end 

of  the  time  interval,  respectively.  Hie  total  work  done  on  the  system 
at  time  t added  to  the  initial  kinetic  energy  of  the  first  particle 
was  then  compared  with  the  total  internal  energy  of  the  lattice  and  the 
two  results  were  found  to  agree. 
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Finally,  we  repeated  a large  number  of  our  calculations  varying  the 
step  size  At  and  consistently  found  that  similar  results  were  produced. 
It  was  found,  however,  that  as  increased,  it  was  necessary  to 

decrease  the  step  size,  At,  from  values  acceptable  for  the  harmonic 
case . 


5.  PROPAGATION  IN  THE  INITIALLY 
QUIESCENT,  ANHARMONIC  LATTICE 

In  this  section,  we  shall  be  primarily  concerned  with  the  numerical 

solution  of  Eq.  (4.9)  for  various  values  of  the  nonlinearity  parameter, 

A . Prior  to  excitation  by  the  shock,  each  of  the  atoms  in  the  lattice 
m ^ 

is  at  rest  in  its  equilibrium  position.  In  the  long-time  limit,  we 
observe  well-defined,  stable  pulses  (.solitary  waves)  of  constant  ampli- 
tude propagating  in  the  vicinity  of  the  shock  front.  These  solitary 
waves  are  rather  unusual  physical  entities  and  we  digress  in  Sec. 

S.J  to  discuss  their  characteristics  and  physical  significance.  In 
Sec.  5.3,  we  present  an  analytical  approximation  for  the  shape  and 
amplitude  of  the  solitary  waves  under  certain  limiting  conditions  and 
compare  the  results  with  our  numerical  solutions.  Finally,  we  conclude 
this  section  by  briefly  considering  the  results  of  our  calculations  for 
shock  propagation  in  a lattice  whose  atoms  interact  via  the  hard-sphere 

potential  of  Northcote  and  Potts  discussed  in  Sec.  2. 

5.1  The  .Morse  Interaction 


We  have  solved  Eq . (4.9)  for  various  values  of  the  nonlinearity 
parameter,  A^,  which  is  a representation  of  both  the  strength  of  the 

compression  and  the  anharmonicity  of  the  lattice.  [See  Eq.  (4.10).] 
The  results  of  the  calculation  can  be  understood  most  easily  by  com- 
paring the  velocity-time  trajectories  for  various  particles  in  the 
lattice  and  indicating  what  conclusions  can  be  drawn  about  the  shock 
profile  in  general. 

Velocity-time  trajectories  for  several  particles  in  the  lattice 
arc  shown  in  Fig.  5 for  the  case  in  which  the  nonlinearity  parameter, 
A^,  was  equal  to  0.2.  In  order  to  plot  these  trajectories  on  the  same 

time  scale,  we  have,  in  each  case,  readjusted  the  time  axis  so  that  at 
t = 0 the  particle  under  consideration  first  feels  the  effect  of  the 
shock  Front.  Ilje  actual  time  at  which  the  front  reaches  the  particle, 
denoted  by  t^^,  is  indicated  in  the  figure. 

For  a particle  close  to  the  driven  end  of  the  chain  (jiarticle  5 
in  the  figure),  the  velocity  variation  is  similar  to  that  found  for 
the  harmonic  lattice.  (See  Fig.  4,  but  note  the  change  in  scale.) 


Ill 


Figure  5.  Vcloci ty-time  trajectories  for  the  initially  quiescent.  Morse-potential  lattice  for 


Upon  cxcitut ion  by  tht-  shock  wave,  the  particle  velocity  reaches  an 
aitplituile  ot'  about  1.4J  anti  subset|uent ly  oscillates  about  unity,  the 
compression  velocity,  mth  ilecreasing  amplitude.  By  the  time  the 

shock  front  has  reached  the  25^^  particle,  the  magnitude  of  the 
oscillations  has  increased.  Meanwhile,  the  pulses  have  become  narrower 
and  more  well  defined.  Physically,  the  early-time  development  of  the 
os,  1 1 l.it  ions  IS  apparently  governed  by  dispersion,  but  as  the  shock 
front  propagates  farther  into  the  lattice,  the  nonlinear  effects 
become  increasingly  import.int.  These  effects,  of  course,  tend  to 
steepen  the  pulse.  It  is  interesting  that  the  two  effects  become 

ir.i-Mtant  in  the  reverse  order  from  that  observed  by  Zabusky  and  Kruskal 
in  their  study  of  oscillations  in  a plasma. 

By  the  time  the  tront  has  reached  the  100^^  particle  the  shape  of 
the  [mlses  has  become  .still  more  well  defined  and  the  amplitude  appears 
to  be  apjiroach i tig  a constant  ma.ximum  value  of  about  2.  If  we  view  a 
p.  •ticiitar  point  in  the  lattice,  the  pulses  in  the  vicinity  of  the 
shock  front  appeir  to  propagate  into  the  lattice  from  the  compressed 
end.  Since  it  is  found  that  the  velocity  with  which  they  propagate 
decreases  with  decreasing  amplitude,  they  tend  to  sp^^ad  apart  as  they 
prop.ig.ite.  When  the  shock  front  has  reached  the  700  particle,  we 
find  that  the  leading  pulses  have  evolved  into  a senuencc  of  e.xtremely 
well-defined  e.xci  tat  ions  (solitary  waves)  and  have  indeed  reached  a 
constant  amplitude  of  approximately  Z.  The  pulse^l^are  more  distantly 
spaced  in  time  than  when  the  front  was  at  the  100  particle,  owing 
to  the  spreading  effect  discussed  previously. 

Therefore,  for  = 0.2,  the  velocity  trajectory  for  particles 

with  N > 200  can  be  described  as  follows:  The  trajector>'  initiall}’ 

varies  along  a sequence  of  solitary  waves  of  essentially  constant 
sha[ie  and  amplitude,  the  number  of  these  waves  increasing  slowly  with 
increasing  particle  number.  I'he  velocity  then  undergoes  damped 
oscillations  about  the  value  unity.  Iventually  these  oscillations  bi'- 
come  negligible  .ind  the  paitiLle,  for  .ill  practical  purposes,  then  moves 
at  a constant  velocity  equ.il  to  the  compression  velocity  for  all  later 

20 

times.  This  bt'h.iviur  is  simil.ir  to  th.it  observed  bv  Tasi“  in  his  studi 
of  a weakly  nonlinear  lattice  described  by  a cubic  interatomic  potent i,i 

Figure  b shows  .i  simil.ir  sequence  of  trajectories  for  the  cast'  in 
which  the  nonlinearity  p.irameter,  A^,  w.is  equal  to  unity.  Qu.i  1 i t .i  1 1 ve  1 \- 

.^4.  .N..J.  Zabusky  and  '■1.0.  kruskal,  "Interaction  of  Solitons  in  a 

(lol  1 i s i on  1 ess  Plasma  and  the  Recurrence  of  Initi.il  .States", 

I'hys.  Rev.  Tetters  l.'i,  240  ( 1 0(>5 ) . 
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Figure  6.  Velocity-time  trajectories  for  the  initially  quiescent,  Morse-potential  lattice  for 
the  case  A =1.0. 
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the  siime  description  applies  near  the  front  except  that  the  pulses 
are  narrower,  owing  to  the  increased  nonlinearity.  A rather  remarkable 
result  of  the  calculations  is  that  the  solitary  waves  evolve  to  the 
same  maximum  amplitude  for  any  nonzero  value  of  A^;  the  time  necessary 

for  the  solitary  waves  to  develop,  however,  increases  with  decreasing 

values  of  A . 

m 

For  A^  = 1.0,  we  again  find  that  the  solitary-wave  train  is 

followed  by  a region  where  the  velocity  oscillates  with  decreasing 
amplitude  about  unity.  In  this  case,  however,  the  oscillations  do  not 

decay  entirely,  but  appear  to  approach  a constant  amplitude.  This 

behavior  was  noted  by  Strenzwilk^^  in  calculations  performed  with  the 
model  discussed  in  this  report.  We  illustrate  this  feature  by  plotting 
in  Fig.  7 the  maximum  velocity  for  particles  behind  the  front  for  both 
values  of  the  nonlinearity  parameter  at  a time  when  the  front  is  approxi- 
mately at  the  480*^^  particle.  The  minimum  velocity  for  each  case  can 
he  obtained  by  reflecting  the  corresponding  curve  about  the  value  unity. 
It  is  apparent  from  the  figure  that,  for  A^  = 1.0,  the  velocity  of 

particles  near  the  end  of  the  compressed  region  oscillates  between 
the  values  1.54  and  0.46.  We  observed  no  tendency  for  these 
oscillations  to  decay  during  the  course  of  the  computation.  This 
result  suggests  that  the  oscillations  far  behind  the  front  approach 
a constant  amplitude  which  decreases  with  decreasing  values  of  the 
nonlinearity  parameter.  For  A^  = 0.2,  this  amplitude  is  so  small 

that  the  particles  are  essentially  moving  at  a constant  velocity 
equal  to  the  compression  velocity. 

I'he  question  arises  as  to  whether  these  oscillations  might  lead 
to  a temperature  rise  behind  the  front,  [hiring  the  calculation,  there 
appeared  to  be  no  tendency  towards  randomization  of  the  velocity;  in 
fact , tive  velocity  of  each  particle  continued  to  vary  along  well-defined 
oscillations.  The  absence  of  thermalization  is  illustrated  in  Fig.  8 
which  shows  the  velocity  distribution  function,  f,  for  the  first  250 
particles  in  the  lattice  at  two  slightly  different  times  for  A^  = 1.0. 

I'he  distribution  function  represents  the  number  of  particles  with 
velocities  in  a particular  interval,  normalized  to  the  total  number  of 
particles.  For  the  system  to  be  in  thermal  equilibrium,  the  distri- 
butio;i  should  be  constant  in  time  and  the  profile  should  correspond  to  a 
Maxwellian  distribution  about  unity  corresponding  to  the  equilibrium 
temperature.  It  is  apparent  from  the  figure  that  the  distribution 


.55.  n.F.  .Strenzwilk,  "Lffect  of  Different  Initial  Accelerations  on  the 
■Subsequent  Shock  Profile  in  One-Dimensional  l.attices",  BRl,  Report 
fto  be  publishedl . 
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Figure  Maximun  particle  velocity  behind 


function  for  this  case  is  not  Maxwellian.  In  fact,  the  shape  of  the 
distribution  curve  changes  dramatically  with  time,  depending  on  where 
in  the  oscillation  cycle  most  particles  happen  to  be  when  the  distri- 
bution is  sampled.  We  conclude,  therefore,  that  shock  compression  in 
the  initially  quiescent  lattice  does  not  lead  to  a temperature  rise 
behind  the  front. 

The  conclusions  which  can  be  drawn  from  our  calculations  as  well 

as  preceding  studies^^  of  shock  propagation  in  an  anharmonic, 
quiescent,  one-dimensional  lattice  can  be  summarized  as  follows:  (1) 

The  shock  speed  is,  on  the  average,  constant  and  increases  with 
increasing  A^.  The  calculation  is  quite  simple  from  the  computer  data 

and  will  not  be  carried  out  here.  (2)  The  early  development  of  the 
shock  profile  is  governed  by  dispersion  which  leads  to  an  oscillatory 
velocity  profile.  As  these  oscillations  grow  in  amplitude,  resulting 
in  larger  relative  displacements,  nonlinear  effects  become  important. 
F.ventually  a point  is  reached  at  which  the  nonlinear  effects,  which 
tend  to  steepen  the  pulse,  exactly  balance  the  dispersion,  which 
tends  to  broaden  it.  Solitary  waves  propagate  in  the  vicinity  of  the 
front  at  this  time.  (S)  The  effects  of  nonlinearity  become  more 
pronounced  with  increasing  A^  for  a given  distance  into  the  lattice, 

and  with  increasing  distance  into  the  lattice  for  a given  A^.  These 

effects  are  made  evident  in  two  ways.  The  solitary  waves  develop  more 
quickly  in  more  nonlinear  lattices,  and,  for  a given  A^,  the  number 

of  solitary  waves  increases  as  the  shock  propagates  farther  into  the 
lattice.  (4)  The  energy  deposited  into  the  lattice  by  the  shock  front 
appears  either  as  potential  energy  due  to  compression  of  the  lattice 
or  as  ordered  translational  energy  of  the  particles  behind  the  front. 
There  is  no  thermal i zation  of  the  energy  so  that  the  compressed  lattice 
remains  at  zero  temperature  after  the  shock  has  passed.  (5)  The  shock 
profile  is  not  steady  in  time  because  of  the  spreading  effect  of  jHilses 
of  different  amplitude  discussed  previously.  The  last  two  results  are 
in  obvious  contradiction  to  the  usual  assumptions  in  continuum  treat- 
ments of  shock  propagation. 

5 . 2 Characteristics  and  Physical  Significance  of  Solitary  Waves 

.Solitary  waves,  such  .as  the  ones  discussed  in  the  last  subsection, 
have  been  observed  to  propagate  in  other  nonlinear,  dispersive  media, 
rhey  represent  particular  solutions  to  nonlinear,  differential  equa- 
tions which  describe  wave  propagation  in  certain  continuous  media  and 

in  the  so-called  "Toda"  discrete  lattice'^^’’'^  . Two  excellent  review 


So.  M.  I'oda  , "Vibration  of  a Chain  With  Nonlinear  1 nt  er.ac  t ion" , .1. 
I’hys . Sue.  .lapan  4.S1  {19b7). 

S7,  M.  Toda,  "Wave  Propagation  in  Anharmonic  l.attices",  .1.  I’hys.  Soc . 
Japan  .SOI  (19b7). 
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articles  ’ have  recently  been  written  which  discuss  the  properties 
of  solitary  waves  and  the  sort  of  systems  to  which  they  apply.  The 
article  by  Toda  emphasizes  in  particular  propagation  in  the  nonlinear 
lattice. 

Until  recently,  solitary  waves  were  considered  to  be  only  of 
academic  interest.  Since  they  represented  only  particular  solutions 
to  differential  equations,  it  was  suspected  that  they  would  exist 
only  under  a rather  special  set  of  initial  conditions.  Consequently, 
it  was  felt  that  their  significance  to  an  arbitrary  initial-value 
problem  was  minimal.  At  any  rate,  it  was  believed  that,  upon  collision, 
two  solitary  waves  would  scatter  irreversibly,  so  that  their  presence 
would,  at  most,  be  a transient  effect. 

The  last  assumption  was  disproved  by  Zabusky  and  kruskal  in 
their  numerical  solution  of  the  Korteweg-de  Vries  equation  describing 
plasma  waves.  Solitary  waves  of  various  velocities  were  allowed  to 
collide  and,  after  collision,  were  observed  to  maintain  their  original 
shapes.  The  term  "solitons"  was  coined  by  Zabusky  and  Kruskal  to 
refer  to  solitary  waves  which  remained  stable  upon  collision  with  one 
another.  The  most  significant,  and  rather  startling,  conclusion  from 
these  investigations  is  that,  for  a system  in  which  solitons  propagate, 
one  can  expect  "randomization"  of  the  energy  or  thermal  equilibrium 
to  occur  only  after  extremely  long  times,  if  ever. 

Zabusky  and  Kruskal  suggested  that  the  existence  of  solitons  could 
explain  the  paradoxical  results  obtained  in  the  famous  i'ermi  , I'ast.i, 

IJlam  computer  study^  in  the  1950's.  These  investigators  had  excited  a 
single  normal  mode  in  a weakly  anharmonic  lattice  and  had  observed  th.it  . 
rather  than  becoming  equipartitioned  among  the  various  normal  modes, 
the  energy  flowed  periodically  among  only  certain  of  the  modes.  In 
addition,  after  some  recurrence  time,  the  system  returned  to  its  initial 

40 

state,  apparently  disproving  the  ergodic  hypothesis  of  statistic.il 

4 1 - 4 .3 

mechanics.  Various  explanations  ‘ of  the  results  have  been  offered. 


,38.  A.C.  Scott,  F.Y.F.  Chu  and  D.W.  Mcl.aughlin,  "The  Soiiton:  A \e . 

Concept  in  Applied  Science",  Proc . IFTili  61,  144.3  (lO'Si. 

.39.  M.  Toda,  "Studies  in  a Nonlinear  Lattice".  Phvsics  Reports 
18C,  1 (1975). 

40.  R.C.  To  1 man.  The  Principles  of  Statistical  Mechanics  (Oxford 
University  Press,  New  York,  19.38),  Chap.  .3. 

41.  ,J.  Ford,  "Fquipart i t ion  of  (-nergy  for  Nonlinear  Systems",  ,1. 
Math.  Phys.  2,  .387  (1961). 

42.  .1.  Ford  and  ,1.  Waters,  "Computer  Studies  of  Fnergy  Sharing  .ii>  i 
Frgodicity  for  Nonlinear  Oscillator  Systems",  ,1.  Math.  I’hvs. 

4,  129.3  (196.3). 

4.3.  r;.A.  .Jackson,  "Nonlinear  Coupled  Oscillators.  I.  Perturbation 

Theory;  Lrgodic  Problem",  .1.  Math.  Phys.  4,  551  (19(i.3). 
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but  the  suggestion  of  Zabusky  and  Kruskal  seems  to  be  the  most  con- 
clusive. Since  their  study,  there  has  been  renewed  interest  in  the 
study  of  solitary-wave  behavior,  as  can  be  witnessed  from  the  review 
articles  cited  previously. 


S . Analytic  Approx  imation  for  Solitary  haves  in  the  Continuum  Limit 


A comparison  of  the  solitary-wave  profile  for  A^  = 0.2  (Fig.  51 

with  that  for  A =1.0  {Fig.  6)  indicates  that  the  solitary  wave 
becomes  broader  as  the  nonlinearity  parameter  decreases.  This  trend 
suggests  that  for  sufficiently  small  values  of  A^  it  might  be  possible 

to  derive  an  analytic  approximation  to  the  solitary-wave  profile  from 
the  continuum  limit  of  the  equations  of  motion.  Therefore,  we  nov 
consider  the  case  in  which  the  nonlinearity  parameter  is  small  and 
take  the  continuum  limit  of  Hu.  (4.91.  For  small  values  of  A , the 

exponential  terms  in  Hq.  (4.9)  can  be  expanded,  and  if  we  retain  onl> 

the  lowest-order  nonvanishing  terms  in  A , we  obtain 
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(5.  1 ) 


The  right-hand  side  of  liq.  (.5.1)  is  just  the  force  derived  from  the 
cubic  interatomic  potential  used  by  Tasi^"^  and  by  Manvi  et  al.^^’ 


We  now  assume  that  can  be  related  to  Sj^  by  a Taylor  expansion 

about  the  lattice  constant,  equal  to  zero,  i.e.,  a continuum  ajqiroxi- 

mation.  Letting  Sj,  = S,  we  have 
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.Substituting  Hq.  (5.2)  into  H.q.  (5.1)  and  retaining  terms  through  the 
fourth  derivative  produce  the  result 
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59 

liquation  (5.. 5)  is  the  so-called  Zabusky  equation  " . If  we  neglect  the 
last  two  terms,  we  get  the  ordinary  linear  wave  equation.  (Recall  that 
the  sound  speed  in  the  harmonic  lattice  was  one-half  in  our  nondimen- 
sional  units.)  The  last  term  on  the  right-hand  side  represents  the 
lowest-order  nonvanishing  term  in  the  nonlinearity,  whereas  the  fourtli- 
dcrivative  term  must  he  kept  to  account  for  dispersion. 
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Equation  (5.3)  is  obviously  difficult  to  solve  in  general,  but  since 
we  are  only  interested  in  determining  what  well-defined  pulses  the  lattice 
can  support,  we  assume  a traveling-wave  solution  of  the  form 


S = S(k  - Ct)  = S(y) 


(5.4) 


where  C is  the  solitary-wave  propagation  velocity.  Substitution  of 
Eq.  (5.4)  into  Eq.  (5.3)  and  noting  that 


yield 
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where  the  primes  denote  differentiation  with  respect  to  y.  Equation 
(5.7)  can  be  integrated  once  and,  applying  the  boundary  conditions  that 
V and  its  derivatives  vanish  at  ± we  obtain 


12C(4C"  - 1)V  = CV''  + 18  A 


(5.8) 


Multiplying  Eq.  (5.8)  by  the  integrating  factor,  V'',  and  integrating 
once  more  yield 


120(40"^  -1)V^  = C(VM"  + 12  A , 


whence 

V' 


= p’(4C“  - 1)V“  - 12  A^  V^/C  j . 

Equation  (5.10,  can  be  solved  for  y to  produce 


(5.9) 


(5.10) 


y = 


f dV 

J j^l2(4(r  - 1)V“  - 12  A^  V^/C  j 


(5.11) 


The  integration  in  Eq.  (5.11)  is  elementary  and  produces  the  result 

V = sech^  ^ 3(4C“  - 1)  y | . (.S.12) 

m 


SO 


liquation  (5.12)  should  represent  an  approximation  to  the  amplitude 
and  shape  of  the  solitary  waves  for  sufficiently  small  values  of  A 

provided  the  continuum  approximation  is  valid.  To  check  the  results, 
we  obtained  from  our  computer  data  values  of  the  propagation  velocity  C 
for  different  values  of  A^,  substituted  them  into  Eq.  (5.12),  and  compared 

the  results  with  the  results  of  our  numerical  calculations.  The  results 
are  indicated  in  Table  II  in  which  we  show  the  value  of  A^;  the  propa- 
gation velocity,  C;  the  soliton  amplitude.  A;  and  the  full  width  of 
the  soliton  at  half  ma.ximum,  A.  The  last  two  columns  show  the 
corresponding  values  as  calculated  from  Eq.  (5.12)  using  values  for  C 
contained  in  column  2.  The  results  are  in  fair  agreement  (within  about 
15%)  for  the  smallest  value  of  A^,  but  the  difference  between  computed 

and  analytical  results  increases  rapidly  with  increasing  A^.  Figure  9 

shows  a comparison  of  the  numerical  and  analytical  wave  profiles  for 

A =0.1. 
m 


TABLE  II.  ANALYTICAL  AND  COMPUTATIONAL  VALUES  OF  SOLITARY-WAVE 

AMPLITUDES  AND  WIDTHS.  VELOCITIES  ARE  EXPRESSED  IN  NUMBER 
OF  PARTICLES  PER  UNIT  OF  NONDIMENSIONAL  TIME. 


Computational  Analytical 


A 

m 

C 

A 

A 

A 

A 

0.1 

0.59 

2.0 

3.1 

2.3 

2.8 

0.2 

0.68 

2.0 

2.3 

2.9 

1.6 

1.0 

1.19 

2.0 

1.0 

5.6 

0.20 

5 . 4 The  Hard-Sp here  Mode 1 of  Northcote  and  Potts 

■>S 

Northcote  and  I’otts  ' have  suggested  a potential  in  which  the  atoms 
are  assumed  to  have  a finite  diameter,  d,  and  interact  with  one  another 
via  harmonic  forces.  Between  collisions,  the  motion  is  equivalent  to 
that  of  a harmonic  lattice,  but,  upon  collision,  the  atoms  exchange 
momenta.  This  extremely  nonlinear  potential  was  used  by  Northcote  and 
Potts  to  repeat  the  Fermi,  Pasta,  Ulam  problem  in  an  effort  to  determine 
whether  the  strong  nonlinearity  would  cause  the  system  to  approach 
thermal  equilibrium.  Their  results  suggested  that  the  system  did 
equilibrate  after  excitation  of  a single  mode.  That  thermal i zat ion  of 
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ot  the  energy  occurred  (perhaps)  led  Montroll  to  speculate  that  this 
potential  did  not  possess  solitary-wave  solutions. 

We  have  modified  our  computer  code  to  treat  a hard-sphere  lattice 
subjected  to  shock  compression.  Our  interest  was  to  determine  both 
whether  equilibration  of  the  energy  occurred  behind  the  front  and 
whether  the  energy  could  be  propagated  in  the  form  of  solitary  waves. 

For  purposes  of  calculation,  we  assumed  a nondimensional  diameter  (real 
diameter  normalized  by  u/u^)  of  0.75  and  a nondimensional  lattice  spacing. 


R 


^o'^o 


u 


(5.13) 


of  3.0.  Again,  we  assumed  that  the  atoms  of  the  lattice  were  at  rest 
prior  to  excitation  by  the  shock. 

Figure  10  shows  the  velocity-time  trajectories  of  the  25^^  and  50^^ 
particles  in  the  lattice  subsequent  to  excitation  by  the  shock.  In  the 
vicinity  of  the  front,  we  find  again  a series  of  well-defined  pulses 
having  constant  amplitude  and  width.  The  amplitude  is  approximately  1.6. 
It  should  be  observed  that  well-developed  solitary  waves  are  propagating 

in  the  vicinity  of  the  front  as  early  as  the  particle.  The  increased 

nonlinearity  causes  the  solitary  waves  to  form  more  rapidly  than  in  the 
Morse-potential  cases  we  studied. 


From  Fig.  10,  it  may  be  noted  that,  at  any  given  particle,  there 
appears  to  be  a tendency  for  the  solitary  waves  to  become  more  distant I>- 
spaced  in  time  with  increasing  time.  It  is  interesting  that  this  is  the 
opposite  effect  from  that  observed  in  the  case  of  the  Morse  potential. 

The  spreading  can  be  seen  more  easily  in  Fig.  11  in  which  is  plotted  the 
velocity  as  a function  of  particle  number  at  three  different  times.  Tlio 
figure  indicates  that  the  length  of  "dead  spaces"  or  regions  of  the  lattic 
in  which  the  particle  velocity  is  essentially  zero  increases  with  both 
distance  from  the  front  and  with  time.  In  addition,  there  appears  to  be 
no  gradual  decay  of  the  pulse  amplitude  as  we  approach  the  region  of 
uniform  motion  behind  the  front,  but  rather  there  is  a sudden  transition 
between  this  region  and  the  region  of  solitary-wave  propagation. 

The  reasons  for  the  differences  in  detail  between  the  hard-sphere 
case  and  the  Morse-potential  case  are  presently  not  clear  to  us.  It  is 
nevertheless  interesting  that  both  predict  a nonsteady  profile  and  no 
temperature  rise  behind  the  front  as  can  be  observed  from  Fig.  11.  In 
addition,  the  hard-sphere  potential  docs  appear  to  be  able  to  support 
the  propagation  of  solitar)’  waves  at  least  under  this  rather  ordered 
set  of  initial  conditions. 


44.  F;.  Montroll,  in  "Proceedings  of  the  bxplosives  (ihcmical  Reactions 

Seminar",  ARO-I)  Report  70-1,  Durham,  NC,  1968,  p.  145. 


6.  PROPAGATION  IN  Tllli  ANTURMONIC  lAn'ICP 
AT  NONZPRO  INITIAL  TEMPERATURn 


r 


The  calculations  of  the  preceding  section  represent  somewhat  of  a 
compromise  to  physical  reality  because,  as  we  have  emphasized,  it  was 
assumed  that  the  particles  were  initially  at  rest  in  their  equilibrium 
I)ositions.  Consequently,  the  effect  of  nonzero  ambient  temperature 
upon  the  shock  profile  has  not  yet  been  examined.  In  this  section, 
we  discuss  calculations  performed  for  the  case  in  which  the  particles 
in  the  lattice  were  allowed  to  undergo  some  initial  thermal  oscillations 
prior  to  being  excited  by  the  shock.  Our  intention  is  to  determine  which, 
if  any,  of  the  anomalous  results  of  the  preceding  section  (nonsteady 
profile,  absence  of  equilibration,  propagation  of  solitary  waves)  are  a 
manifestation  of  the  rather  ordered  set  of  initial  conditions  used  in 
that  calculation. 

In  Sec.  6.1,  we  discuss  the  procedure  for  doing  the  calculations 
and  in  Sec.  6.2  the  method  for  preparing  a one-dimensional  lattice  iii 
thenrial  equilibrium.  In  Sec.  6.,T,  various  velocity-time  trajectories  arc 
discussed  as  before  and  we  observe  high-amplitude  solitary  waves  propa- 
gating amid  a thermal  background.  The  solitary  waves  are  isolated  from 
the  thermal  background  and  compared  with  those  obtained  in  the  previous 
section.  In  Sec.  6.4,  we  observe  a "collision"  or  nonlinear  interaction 
between  two  solitary  waves  and  discuss  the  consequence  of  their  stability. 

In  Sec.  6.5,  we  investigate  the  possibility  of  thermal  equilibrium 
behind  the  shock  front  and  calculate  the  thermodynamic  variables  of  interest . 
Finally,  in  Sec.  6.6,  we  examine  the  extent  to  which  the  steady-state 
conservation  equations  are  satisfied. 

6.1  Niethod  for  Pejformr^  Calculat  ions 

Since  the  shock  front  travels  at  a finite  speed  into  the  lattice 
and  since  our  major  interest  is  in  the  shock  profile,  it  is  clearly 
unnecessary  to  monitor  the  particles  which  lie  well  ahead  of  the  front. 

In  the  preceding  section,  this  posed  no  probiom  because  the  particles 
ahead  of  the  front  were  at  rest  in  their  equilibrium  positions.  Con- 
sec(uently,  we  simply  provided  our  computer  pi'ogram  with  a test  which 
caused  the  calculation  to  include  only  particles  whose  velocity  and 
displacement  were  greater  than  some  small  number. 

A similar  test  is  no  longer  possible  in  the  present  calculation 
because  the  particles  ahead  of  the  front  arc  now  in  motion  prior  to 
being  excited  by  the  front.  In  order  to  avoid  monitoring  all  the 
particles  ahead  of  the  front,  it  is  clearly  desirable  to  add  the 
particles  in  segments  only  as  they  are  needed.  One  might  still  exjiect , 
however,  that  it  would  be  difficult  to  detect  the  location  of  the  shock 
front  in  the  lattice  particularly  for  we.ik  shocks. 


An  extremely  clever  solution  to  the  problem  was  suggested  to  us 
45 

by  Tsai  . We  simulate  the  semi-infinite  lattice  at  t = 0 by  a series 
of  identical  segments,  each  containing  n particles  in  thermal  equili- 
brium, as  shown  in  fig.  12.  In  other  words,  the  velocities  and  dis- 
placements of  the  particles  in  the  chain  are  initially  related  by  the 
periodic  conditions 


V CO)  = V.  (0) 
m jn+m 


S (0)  = S.  (0) 
m^  jn+m^ 


1 ~> 


m = 1 , 2 , . . .n 


The  first  particle  is  then  subjected  to  steady  compression  at  a 
nondimens ional  velocity  of  unity  as  before  and,  as  a result,  a shock 
wave  proceeds  through  the  chain. 


Consider  the  particle  assumed  to  be  located  in  the  first  segment. 
Before  the  shock  front  reaches  this  particle,  its  motion  is  identical  to 

the  j+n^^'  particle  because  these  atoms  are  acted  upon  by  identical  forces 
and  had  identical  initial  conditions.  Working  our  way  backwards  from 

the  particle,  we  locate  the  first  particle  whose  displacement  and 

velocity  differ  appreciably  from  those  of  the  s>'mmetric  particle  in  the 
second  segment  and,  thereby,  determine  exactly  the  locat  ioj^^^of  the  front 
at  any  time  i.  Furthermore,  until  the  shock  reaches  the  n particle,  it 
is  necessary  to  solve  the  equations  of  motion  for  only  the  first  two 
segments  since  the  particles  in  all  the  remaining  segments  will  have 
trajectories  identical  to  the  corresponding  particles  in  the  second^j^ 
segment.  (In  order  to  calculate  the  force  from  the  right  on  the  2n 

St 

particle,  the  displacement  of  the  2n+l  particle  is  required.  However, 
we  do  not  have  to  monitor  this  particle  since  its  displacement  is 

equivalent  to  that  of  the  n+l^*”  particle.)  When  the  shock  front  nears 

the  n'^*^  particle,  the  particles  in  the  third  segment  are  included  in  the 
computation.  The  shock  front  is  located  in  the  same  manner  as  before, 
and  the  process  repeated  as  often  as  necessary  to  complete  the  calculation. 
It  is  therefore  necessary  to  monitor  at  most  2n  particles  ahead  of  the 
shock  front  at  any  time,  the  last  n representing  the  equilibrium  condition- 
of  the  chain  ahead  of  the  shock.  For  purposes  of  calculation,  n was 
taken  to  be  100. 


() . 2 Preparation  of  Initial  S e gment  in  lltermal  F.qu  i 1 ihr  ium 


We  begin  this  subsection  by  considering  the  properties 
dimensional  chain  of  atoms  in  thermal  equilibrium.  We  then 
we  may  artificially  prepare  tlie  lattice  so  these  properties 


of  a one- 
discuss  how 
are  sat  i sf  ii\i . 
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Figure  12.  Construction  of  a semi-infinite  chain  from  a sequence  of  initially  identical  segments. 


Consider  a chain  of  n particles  obeying  cyclic  boundary  conditions 

such  that  the  first  and  n*^^  atoms  are  connected.  The  chain  is  assumed  to 
be  in  thermal  equilibrium  at  temperature  T.  From  the  definition  of 
equilibrium,  the  probability  that  the  system  occupies  a particular 
differential  volume  in  phase  space  is  given  by  the  distribution  function 

r 2 

/ mv  . \ 

F(Vj,...v^,  Xj,..,  Xj^)dv^.  . .dv^dxj. . .dx^  = ^ exp  j /l^'I 

dv, . . .dv  dx, . . .dx  (6.11 

Inin 


where  2 is  the  partition  function  which  normalizes  the  distribution  to 
unity  and  4>  is  the  total  potential  energy  of  the  lattice.  The  probability 
that  particle  j has  velocity  between  v^  and  v^+dv^  regardless  of  its 

displacement  and  regardless  of  the  velocities  and  displacements  of  the 
remaining  particles  of  the  lattice  is  obtained  by  integrating  Hq.  (6.1) 

over  all  displacements  and  all  velocities  except  the  i’"*'.  If  we  assume 
a velocity-independent  potential,  the  result  is 


F(v^)dVj  = ^ exp  - mv"j/(2kT) j dv^ 


where 


exp  [-  mv‘'^/(2kr)]  dVj 


liquation  (6.2)  indicates  that  the  velocities  of  all  particles  are 
distributed  according  to  a Maxwellian  distribution  at  temperature  T, 
regardless  of  the  form  of  the  potential  interaction  as  long  as  it  is 
velocity  independent.  Consequently,  the  mean  (ensemble  average)  square 
velocity  of  any  particular  particle  is  given  by 


< v“  > 


00 

J"  v“  exp  -mv“/ (2kT)  j 


dv  = kT/m 


and  the  mean  thermal  velocity  by 


"t  = y < V-  > 


w 
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If  the  lattice  discussed  above  were  harmonic,  it  is  well  known  that 
the  total  energy  contained  in  the  lattice  is  given  by 

E = nkr  . (b-b) 

On  the  average,  half  of  the  energy  is  potential  and  half  kinetic.  If  an 
initial  energy  given  by  bq . (tt.b)  is  put  into  an  anharmonic  lattice  and 
the  lattice  allowed  to  come  to  thermal  equilibrium,  it  will  equilibrate 
at  some  slightly  different  temperature  owing  to  the  difference  in  the 
harmonic  and  anharmonic  potentials.  If  the  lattice  is  only  slightly 
anharmonic,  however,  the  temperature  T'  should  not  be  appreciably  differ- 
ent from  T and  F.q.  (b.b)  will  be  a reasonable  estimate  of  how  much  energy 
to  initially  put  into  the  lattice  to  correspond  to  an  anharmonic  lattice 
in  equilibrium  at  temperature  T.  Before  proceeding  further,  it  is 
convenient  to  nondimensional ize  the  quantities  of  interest  as  before.  The 
thermal  velocity,  normalized  by  the  compression  veloci ty, then  is 


and,  the  total  energy  in  the  lattice,  normalized  by  the  kinetic  energy 
of  the  first  particle,  is 


e . 2nkT 
'amu"  mu*" 


(b.H) 


ITio  thermal  velocity,  Vj,,  is  the  only  additional  parameter  that  must  be 
specified  in  the  computer  input. 


In  order  to  prepare  our  initial  segment  of  n particles  in  thermal 
equilibrium,  we  generated  a set  of  n-1  random  velocities  o^^ained  from 
a Maxwellian  velocity  distribution.  The  velocity  of  the  n particle 
was  determined  from  the  condition 


(b.y) 


imposed  so  as  to  impart  no  net  momentum  to  the  chain.  bach  particle  in 
the  ch.iiii  was  .issumed  to  be  in  its  equilibrium  position  initially  (no 
potent lal  energy)  .uul  the  velocit ies  were  scaleil  by  a constant  factor  c. 
iTiis  factor  was  chosen  so  that  the  initial  nondimensional  energy  in  the 
chain  was  given  liy  Iq.  (f'.H).  Ihus,  we  required 


n 


.l  = l 

which  determines  the  factor 


M) 


' 


(b, lO) 


Having  obtained  the  scaled  velocities  in  the  manner  discussed 
above,  we  assigned  them  to  each  of  the  n particles  in  the  segment.  The 
lattice  was  then  allowed  to  oscillate  freely  and  we  employed  a number  of 
checks  to  insure  that  it  was  actually  in  a state  of  thermal  equilibrium. 
First,  we  calculated  the  velocity  distribution  function  at  several  times 
and  found  that  the  function  was  constant  in  time  and  essentially  Maxwellian. 
Second,  we  calculated  the  time  average  of  the  velocity  squared  for 
several  of  the  atoms  in  the  lattice,  obtained  by  evaluating  the  expression 


■) 

vr 

1 


V . “ (T)dr 
J 


(6.11) 


where  At  is  a time  long  compared  with  unity.  We  found  that  the  time 
average  agreed  with  the  average  taken  over  all  particles  in  the  system  at 
any  time,  viz., 

< V-  > = ^ L V ' (x)  . (O.i:) 

J 

'ntird,  we  found  that  the  total  kinetic  and  potential  energies  in  the 
lattice  remained  essentially  constant  in  time.  .As  we  pointed  out  before, 
one  would  expect,  for  a harmonic  lattice,  that, the  kinetic  and  potential 
energies  would  be  equal  and  each  given  by  n v j." . For  the  anharmonic  case, 

this  is  no  longer  true  and  we  found  generally  that  the  kinetic  energy  was 
slightly  larger  than  the  potential.  Finally,  we  calculated  the  local 
potential  and  kinetic  energies  by  averaging  over  ;i  few  (say,2.S)  particles 
within  the  segment  and  found  that  these  were  roughly  constant  in  both 
space  and  time.  .Some  deviations  did  occur,  but  these  were  attributed 
to  the  fact  that  relatively  few  particles  were  contained  in  the  average. 

We  should  point  out  that  we  have  studied  equilibration  of  the  chain 
when  the  ends  were  subjected  to  various  boundary  conditions,  but  have 
found  that  only  the  cyclic  condition  produced  reasonable  results.  l-or 
the  case  of  fixed  end  conditions,  we  found  that  the  unequal  forces 
exerted  (<it  any  given  time)  upon  llie  ends  of  the  chain  by  the  fixed 
particles  tended  to  cause  the  center  of  mass  of  the  chain  to  drift 
first  in  one  direction  and  then  in  the  other.  This  resulted  in  rather 
large-scale  displacements  of  the  interior  particles  from  their  equili- 
brium positions.  I'or  the  case  of  free  ends,  we  found  extremel)’  large 
displacements  of  the  particles  particularly  toward  the  ends  of  the 
chain.  It  is  commonly  accepted  tiiat , for  sufficiently  large  systems, 
the  end  conditions  do  not  affect  the  bulk  properties  of  the  system. 
Apparently,  however,  a chain  of  100  particles  is  not  large  enough  to 
neglect  such  effects.  Hie  gener.il  desirability  of  using  periodic 
boundary  coiulitions  in  computer  studies  of  small  systems  has  been 
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discussed  elsewhere^^. 

We  have  considered  in  our  calculations  various  values  of  the  thermal 
velocity  V.^,.  If  we  choose  as  reasonable  values 

T = 300° k 
- ’S 

m = 10  kg 
u = 500  m/s  , 


^ . (b.131 

mu'"  / 

Consequently,  at  room  temperature  we  do  not  expect  the  thermal  velocity 
to  be  negligible  with  respect  to  the  compression  velocity.  Furthermore, 
as  T decreases,  Vj,  only  decreases  as  T,  and  it  is  likely  that  only  at 

very  low  temperatures  indeed  can  one  expect  that  neglecting  the  thermal 
velocity  is  a realistic  assumption.  It  was  this  consideration  wiiich 
prompted  the  calculations  of  this  section. 

t> . 3 Veloc  ity  - I'i r . i j e c tories  and  Prop  a g ;U  ion  of  Solitary  Waves 

In  this  subsection,  we  discuss  the  ve loc i ty-t ime  trajectories 
obtained  when  the  equilibrated  lattice  discussed  above  was  subjected 
to  shock  compression  for  the  case  in  which  = 1.0  and  V.j,  = 0.5. 

Figure  13  shows  the  trajectories  for  three  particles  in  the  lattice, 
namely,  100,  JOO,  and  300.  Again,  the  time  scales  have  been  adjusted 
so  that,  at  time  r = 0,  the  particle  first  feels  the  effect  of  the  shock. 

Referring  to  the  lOo''^  particle  and  to  the  pulses  labelled  A,B,  and  C 
in  the  graph,  we  observe  some  evidence  of  high-amplitude  solitary  waves 
projiagating  ;imid  a thermal  background.  Because  of  the  background, 
however,  the  shapes  of  the  pulses  are  not  so  well  defined  as  in  the  case 
for  which  the  particles  were  initially  at  rest  in  their  equilibrium 
positions.  Recalling  that  propagation  is  occurr in^^toward  the  left,  we 
see  that  by  the  time  the  front  has  reached  the  200^  particle,  pulse  C 
has  begun  to  approach  pulse  B,  whereas  pulse  A has  moved  considerably 
farther  ahead  of  B.  At  particle  300,  pulse  A has  moved  still  farther 
ahead  of  B,  but  it  appears  that  pulses  B and  C ^^e  nearly  the  same 
distance  apart  as  when  the  shock  was  at  the  200  particle. 
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We  are  therefore  led  to  speculate  that  B and  C crossed  while  propagating 
between  the  200^^  and  300^^  particles. 

The  above  observations  are  rather  imprecise.  In  fact,  it  is  not 
really  clear  whether  the  pulses  we  have  observed  are  high-amplitude 
solitary  waves  interacting  with  a thermal  background  or  simply  rather 
extreme,  perhaps  unstable  thermal  oscillations.  In  an  effort  to  resolve 
this  point,  we  redid  the  calculation,  but  instantaneously  stopped  the 

compression  when  the  shock  front  was  at  the  140^^  particle.  The  remaining 

particles  of  the  lattice,  those  beyond  the  140^^,  were  then  placed  at 
rest  in  their  equilibrium  positions.  Essentially,  we  had  a hot,  shock- 
compressed  lattice  directly  next  to  a cold  lattice.  If  the  large-velocity 
oscillations  «ere  indeed  the  result  of  high-amplitude  solitary  waves, 
then  we  should  see  these  waves  propagating  with  superacoustic  speeds 
into  the  cold  lattice  leaving  the  thermal  background  behind.  We  should, 
therefore,  be  able  to  isolate  the  solitary  waves  from  the  thermal 
background . 

Figure  14  illustrates  the  results  of  the  calculation.  The  velocity- 
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time  trajectory  of  the  145  particle  is  plotted  for  the  case  in  which 

the  shock  compression  was  stopped  at  the  140*^^  particle  and  the  remaining 
particles  placed  at  rest  in  their  equilibrium  positions.  For  comparison, 

the  trajectory  of  the  145*^^'  particle  obtained  in  the  original  calculation 
is  shown  in  lower  part  of  the  grapli. 

From  Fig.  14  we  indeed  find  solitary  waves  propagating  out  into  the 
cold  lattice.  Unlike  foe  the  case  studied  in  the  last  section,  however, 
the  pulses  apparently  do  not  evolve  to  the  same  constant  amplitude  and, 
at  the  head  of  the  front,  the  amplitudes  are  higher  than  for  the  case  of 
the  initially  cold  lattice.  Comparison  of  the  location  of  the  peaks  of 
the  p\ilses  in  the  upper  and  lower  parts  of  the  figure  again  emphasizes 
that  the  original  pulsts  which  we  saw  are  solitary  waves  propagating 
upon  a thermal  background.  Consequently,  the  labels  A,B,  and  C, 

correspond  to  the  same  pulses  observed  propagating  between  the  100^^  and 

30()^^part  ic  1 es  in  Fig.  13;  the  label,  1),  refers  to  the  fourth  solitar>’ 
wave  not  discussed  previously. 

We  have  calculated  from  our  computer  data  the  amplitudes,  and 
propagation  velocities  of  waves  A,B,C,  and  U in  the  cold  lattice.  ilie 
results  are  shown  in  Table  III.  We  observe  again  that  the  propagation 
velocity  increases  with  increasing  amplitude.  Consequently,  a spreading 
effect  will  be  introduced  as  before  which  will  prevent  the  profile  from 
approaching  a steady  state.  We  can  also  anticipate  from  the  table  that 
[nilses  B and  C will  eventually  "collide".  A simple  calculation  reveals 

th.1t  the  collision  should  take  place  in  the  vicinity  of  the  280^^^  partule 
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tabu:  111. 


PROP AG A IT ON  VtLOCITY  AND  AMPITTUUK  OF  SOLITARY  WAVES. 
PROPAGATION'  VELOLTTIES  ARE  IN  PARTICLES  PER  UNIT  OF 
NONDIMF.NSIONAL  TIME. 


Sol iton 

Ampl itude 

Propagation  Velocity 

A 

3.2 

1.50 

B 

2.9 

1.42 

C 

3.2 

1.50 

D 

2.7 

1.38 

in  the  lattice  appro.ximately  9b  units  of  time  after  soliton  B reaches  the 

145^^  particle.  Since  the  velocity  of  the  high-amplitude  solitary  waves 
is  not  greatly  influenced  by  the  thermal  background,  our  original  specu- 
lation that  B and  C crossed  prior  to  reaching  the  particle  appears 

correct.  In  the  following  section,  we  examine  the  stability  of  the 
solitary  waves  upon  collision. 

0.4  Soli tary-W ave  Colli  si ons 


Since  we  have  observed  in  the  previous  section  that  the  solitary 
waves  may  collide  at  some  point  in  the  lattice,  it  is  necessary  to 
examine  their  stability  upon  collision.  Such  a collision  is  represented 

in  Fig.  15.  At  the  170^^  particle  in  the  lattice,  we  see  two  solitary 
waves  well  separated  in  time.  Propagation  is  occurring  toward  the  left 
and  the  higher-amplitude,  faster-moving  wave  is  behind  the  slower  one. 

By  the  time  the  propagation  has  reached  the  185^^  particle,  we  find  that 
the  faster-moving  solitary  wave  has  overtaken  the  slower  one  and  a 
nonlinear  interaction  is  occurring.  It  should  be  emphasised  that  the 
interaction  is  nonlinear  and  not  just  a linear  superposition  of  the 
separate  waves  as  might  be  expected  if  the  solitary  waves  represented  the 
solution  to  a linear  differential  e^quation.  When  the  wave  pair  has 

reached  the  200^^  particle,  we  find  that  the  faster  solitary  wave  has 
moved  completely  through  the  slower  one,  but  each  has  maintained  its 
original  shape. 

We  have  monitored  several  of  these  collisions,  and  in  each  case  the 
solitary  waves  emerged  from  the  collision  with  tlieir  original  shapes 
intact.  Consequently,  we  conclude  that  the  solitary  waves  which  we 
observe  propagating  in  the  lattice  are,  in  fact,  solitons.  I'he  presence 
of  these  stable  entities  impedes  the  thennal i zat ion  of  energv  behind  the 
shock  front  and  raises  the  question  as  to  whether  thermal  equilibrium 


at  an  elevated  temperature  is  established  in  the  compressed  lattice, 
iliis  question  will  be  addressed  in  the  following  subsection. 


6. 5 Investigation  of  Thermal  Equilibrium  Behind  the  Front  and 

Calculation  of  the  Thermodynamic  Variables 

We  have  plotted  in  Fig.  16  the  profiles  of  the  potential  and  kinetic 
energies,  averaged  over  100  particles,  for  several  times.  The  single- 
partiole  potential  and  kinetic  energies  were  calculated  from  the  formulas 
indicated  in  Table  I.  The  profiles  exhibit  a structure  characteristic  of 
shock  profiles  in  the  continuum  approximation,  namely,  a transition 
region  joining  two  regions  in  which  the  average  values  of  the  variables 
are  roughly  constant. 

Because  of  the  large  number  of  particles  which  have  been  included 
in  the  average,  it  is  difficult  to  detei’mine  from  the  figure  whether  or 
not  the  transition  region  increases  in  width  as  the  shock  propagates 
farther  into  the  lattice.  However,  it  is  apparent  from  the  discussion 
of  the  preceding  subsection  that  the  transition  region,  which  consists 
of  high-iimpl itude  solitons  propagating  with  different  velocities  on  the 
thermal  background,  must  increase  with  propagation  distance.  The  growth 
of  the  transition  region, owing  to  the  spreading  of  the  sol itons.  should 
become  evident  even  in  the  averaged  profiles  if  the  calculations  are 
extended  to  longer  times.  We  conclude,  therefore,  that  the  shock  pro- 
file is  not  steady  in  time  as  is  generally  assumed  in  continuum  treat- 
ments of  shoi'k  propagation  in  solids. 

The  relatively  constant  energy  densities  located  well  behind  the 
front  lead  us  to  inquire  whether  this  segment  of  the  lattice  is  in 
thermal  equilibrium.  'Hie  question  is  somewhat  difficult  to  answer 
conclusively,  particularly  in  one  dimension,  because  of  the  fairly 
small  number  of  particles  behind  the  front.  If,  for  instance,  we 
calculated  the  velocity  distribution  function,  significant  deviations 
from  a Maxwellian  function  might  bo  anticipated  even  if  the  region 
were  in  f.act  in  equilibrium. 

Nevertheless,  in  an  attempt  to  answer  this  question,  we  located 
from  l ig.  Ih,  at  a particular  time  : , the  appro.\imate  point  behind 
the  front  at  which  the  average  energies  discussed  above  b.'come  const.int  . 
I'he  velocity  distribution  function  for  particles  located  between  this 
point  and  the  driven  end  of  the  lattice  was  then  calculated.  I'he  cal 
dilation  was  subsequently  repeated  at  different  values  of  the  time  .ind 
for  different  numbers  of  particles  contained  within  the  s.imjile.  We 
consistently  obtained  a ne.arly  M.ixwellian  distribution  function  which 
.ippeared  to  cor-espond  to  a thermal  energy  that  w.is  ajiprox  iraat  el  v .i 
f.ictor  of  J.5  higher  th.in  the  ambient  thermal  energy.  A tvpical  dis- 
tribution function  obt.iined  is  shown  in  Fig.  H in  which  the  number  of 
itoms  liaving  velocity  in  .i  given  interval  is  plotted.  The  distribution 
Is  that  obtained  at  time  T=2.B0,  .and  the  sample  consisted  of  the  first 
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potential-  and  kinetic-energy  profiles  for  the  lattice  with  nonzero  initial 
are  at  three  times.  The  solid  and  dashed  lines  represent  the  kinetic  and 
il  energies,  respectively.  The  average  values  are  plotted  at  the  midpoint 

ppropriate  range. 


The  shock  front  at  this  time  was  around 


250  particles  of  the  lattice. 

the  560*'^  particle.  As  can  be  seen  from  Fig.  16,  the  energy  density  is 
nearly  constant  at  this  time  for  this  sample  of  particles.  At  least 
qualitatively,  the  Maxwellian  character  of  the  distribution  function  is 
obvious . 

That  the  distribution  function  remains  constant  in  time  and  corre- 
sponds to  a higher-than-ambient  temperature  (See  Tables  IV  and  V.)  leads 
us  to  conclude,  tentatively,  that  some  equilibration  of  energy  occurs 
behind  the  front.  For  a number  of  reasons,  however,  the  question  has 
not  been  settled  conclusively.  First,  it  is  possible  that  the  times 
at  which  the  distribution  function  was  calculated  were  not  sufficiently 
well  separated  to  discern  time-dependent  behavior  in  the  function.  It  is, 
therefore,  possible  that  at  a much  longer  time  the  function  will  have 
changed  significantly.  Second,  it  is  not  entirely  clear  whether  the 
deviations  from  a theoretical  Maxwellian  distribution  function,  obvious 
in  Fig.  17,  result  from  the  small  number  of  particles  in  the  sample  or 
from  an  actual  state  of  nonequilibrium.  It  is  clear,  however,  that  the 
distribution  function  is  much  more  nearly  Maxwellian  than  those  obtained 
in  Sec.  5 (zero  ambient  temperature),  and  it  seems  unlikely  that  the 
initial  thermal  oscillations  could  have  caused  so  substantial  a change 
in  the  function  without  any  additional  equilibration. 

The  question  of  the  approach  to  thermal  equilibrium  clearly  needs 
further  investigation.  We  will,  however,  defer  such  investigation 
until  the  extension  of  the  calculations  to  throe  dimensions  has  been 
made.  The  larger  number  of  particles  behind  the  front,  in  that  case, 
should  at  least  lead  to  a more  accurate  determination  of  the  distri- 
bution function. 

Assuming  equilibration  far  behind  the  front,  we  have  also  calcu- 
lated other  thermodynamic  variables  both  in  the  compressed  and  uncom- 
pressed, equilibrated  regions  of  the  crystal.  As  we  have  pointed  out 
previously,  the  solution  of  the  equations  of  motion  of  the  particles 
in  the  lattice  is  independent  of  the  equilibrium  lattice  spacing 

However,  in  order  to  calculate  both  the  shock  speed  and  the  densit> 
within  any  region  of  the  crystal,  a value  of  this  parameter  must  be 
specified.  We  chose  a normalized  value,  given  by 


of  6.0  and  found  that  the  shock  speed,  nomialized  to  the  compression  veloc 
ity,  was  essentially  constant  in  time  and  approximately  equal  to  S.(>0.  lliis 
shock  speed  corresponds  to  a Mach  number  of  approximately  5.  Ilu'  pressure 
within  a given  region  of  the  crystal  was  defined  to  be  the  force  exerte^l 
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on  the  particle  by  the  j + averaged  over  the  number  of  particles 
in  the  region  under  consideration.  The  thermal  energy  was  defined 
by  the  expression 

= < (V.  - V)^  > (6.14) 

where  V is  the  average  particle  velocity  in  the  region  and,  here,  the 
notation  < > means  an  average  over  the  appropriate  number  of  particles. 

The  values  of  these  thermodynamic  variables,  as  well  as  the  average 
kinetic  and  potential  energies  discussed  previously,  are  indicated  in 
Tables  IV  and  V.  Table  IV  contains  the  initial  values  or  values  in  the 
uncompressed  region  of  the  lattice.  In  each  case,  these  values  were 
calculated  by  averaging  over  the  last  100  particles  in  the  lattice,  all 
of  which  were  located  ahead  of  the  shock.  In  the  table,  p is  the  nor- 
malized density  in  t!  's  region  of  the  crystal,  or  the  average  number  of 
particles  per  nondimensional  lattice  spacing,  P is  the  average  pressure 
defined  previously,  V is  the  average  particle  velocity,  and  i is  the 
time  at  which  the  calculation  was  made. 


TABLK  IV.  VALUES  OF  THERMODYNAMIC  VARIABLES  IN  UNCOMPRESSED  LATiTCE  Ai 
DIFFERENT  TIMES.  ALL  VALUES  WERE  OBTAINED  BY  AVERAGING  OVER 
THE  lAST  100  PARTICLES  IN  THE  lATTICE. 


T 

V 

kE 

1 

’E 

1 

'T 

P 

P 

:oo 

0 

0.3 

2 

0, 

. 18 

0. 

,32 

0. 10' 

0 . 9(' 

d.SO 

0 

0.2 

8 

0, 

22 

0. 

,28 

0.  lo5 

1.2.' 

I'ABLE  V. 

VAI.UES 

OF 

THERMODYNAMIC 

VARIABLES 

i IN  COMPRESSED 

lATTICI. 

AI 

DIFFER! 

5NT 

ITMES 

AND  FOR 

DU 

•'FERENl 

■ NUMBERS 

OF 

P.\R 

I ICLES 

IN 

THE  SAMI’ LE 

• 

1 

N 

s 

V 

KL 

1 

’E 

E-r 

P 

1’ 

200 

150 

0. 

99 

1.78 

1, 

.01 

0.80 

0. 

197 

5.  O' 

200 

1. 

00 

1.82 

1 , 

,03 

0.82 

0. 

195 

5.72 

250 

150 

0. 

97 

1 .1)5 

1 , 

,07 

0.  T'l 

0. 

197 

0.02 

200 

0. 

93 

1 . 50 

1 . 

, 1 1 

0.  TO 

0. 

I9T 

0.  19 

250 

0. 

92 

1 . 59 

1 . 

,00 

o.:’5 

0. 

194 

5 . SO 

— > 


Results  are  obtained  for  three  different  sample  sizes  by  averaging  over 
the  first  ISO,  first  200,  and  first  250  particles  in  the  lattice.  The 
number  of  particles  in  the  sample  is  indicated  in  the  second  column  of 
the  table  by  the  parameter,  N^. 

Referring  first  to  Table  IV,  we  see  that  the  initial  values  of  the 
various  parameters  are  in  rough  agreement  for  the  two  times  considered. 

It  is  likely  that  more  exact  agreement  could  be  obtained  by  choosing  an 
initial  sample  which  contained  more  than  100  particles  (i.e.,  n larger 
than  100  in  Sec.  6.1),  but  it  is  also  possible  that  we  have  not  achieved 
complete  thermal  equilibrium  in  this  segment.  If  not,  the  deviations 
from  equilibrium  were  felt  to  be  sufficiently  small  for  practical  con- 
siderations . 

It  is  of  particular  interest  to  note  in  Table  V that  the  therm. vl 
energy,  behind  the  front  is  nearly  constant  for  all  samples  and  at 

both  times  considered.  It  is  roughly  a factor  of  2.5  higher  than  the 
ambit- It  thermal  energy  as  we  have  pointed  out  previously.  Since  the 
theritial  energy  is  directly  proportional  to  the  temperature,  the 
temperature  is,  therefore,  also  a factor  of  2.5  higher.  We  observe  an 
increase  in  the  density  of  about  20°  and  an  increase  in  pressure  by 
about  a factor  of  6.  The  average  particle  velocity  in  each  of  the 
samples  is,  of  course,  nearly  equal  to  the  compression  velocity. 

() . 6 The  Conservation  liquations 

If  the  shock  profile  were  steady  in  time,  the  thermodynamic  variaiiles 
in  front  of  and  behind  the  shock  obtained  in  the  last  subsection  would 
satisfy  the  following  conservation  cquation.s: 


p.  (V,  - 11^)  • p,.  (Vf  - Uj.) 


(6. 15a) 


I’i  ^ 4p.  (V.  - iiy)“  = i>^.  + 4p^  (Vj-  - u,,r 


{6 . 15ti  I 


(V,  - P..IV  (V|  - u^)  - (V,,  - 11^1  . 

16. 15c) 


In  these  equations,  denotes  the  shock  speed,  which  is  approximately 

8.(i9  as  discussed  earlier,  and  i and  f denote  initial  values  in  front 
of  the  shock  (Table  IV)  and  final  values  behind  the  shock  (Table  V), 

respectively.  The  ()uantity  li"  is  the  total  energy  measured  in  the 
shock  frame  and  is  n-lated  to  the  potential  and  kinetic  energies  in 
the  lab  fr.ime  through  the  formula 


(6.16) 


= Pt  + Kii  - 2Vl)^  + 

liquations  (6.15)  can  be  derive^',  from  more  general  hydrodynamic  equations’’ 
when  considered  in  conjunction  with  an  equation  of  state,  they  yield  the 

48 

so-called  Hugoniot  conditions 

Although,  as  we  have  stressed  earlier,  the  profile  is  not  steady  in 
time,  it  is  of  interest  to  determine  the  ex.ent  to  which  Eqs.  (6.15)  are 
satisfied.  To  do  so,  w,  i.ccked  the  conse  vation  equations  at  t=250 
using  the  data  from  Table  IV  and  V.  The  final  values  were  taken  to  be 
the  averages  over  the  first  250  particles  in  the  lattice  (last  row, 

I'able  V).  The  equations  were  then  checkt  1 by  directly  substituting  for 
the  variables  in  Eqs.  (o.lS)  and  comparing  the  right-  and  left-hand  sides 
of  the  equations. 

The  results  are  uu'.’cated  in  Table  \1.  Column  1 lists  the  equation 
number,  column  2 the  vali.  ' obtained  for  the  left-hand  side  of  the 
equation,  column  5 the  corresponding  value  for  the  right-hand  side,  and 
column  4 the  percent  deviation.  ITie  last  figure  was  calculated  by  taking 
the  absolute  value  of  the  ratio  of  the  difference  in  the  right-  and  left- 
hand  sides  to  the  value  of  the  left-hand  side.  It  is  rather  interesting 
that  the  largest  deviations  are  around  5%,  suggesting  that  the  nonsteady 
behavior  of  the  profile  affects  the  Hugoniot  conditions  only  slightly. 

A similar  conclusion  has  been  reached  by  Tsai. 


lABl.E  VI.  STEADY-STATE  CUNSIiRVATION  EQUATIONS 
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New  York,  1957),  Chap.  4. 


- 1 


SIJMM/\K\  , CONCLUSIONS,  AND  i UTURL  INVESTIGATIONS 


We  have  carried  out  numerical  and  some  analytic  calculations 
describing  shock  propagation  in  a one-dimensional , discrete  lattice.  A 
number  of  interatomic  potentials  have  been  studied,  but  the  most  exten- 
sive treatment  has  been  given  to  the  case  in  which  the  atoms  interact 
via  a Morse-type  potential,  i'wo  special  cases  of  initial  conditions 
have  been  considered.  In  the  first,  we  assumed  the  atoms  were  initially 
at  rest  in  their  equilibrium  positions  (zero  ambient  temperature)  prior 
to  being  excited  by  the  shock;  in  the  second,  the  lattice  was  initially 
in  thermal  equilibrium  with  an  energy  that  corresponded  roughly  to  room 
temperature. 

ITie  major  conclusions  which  can  be  drawn  from  the  calculations  are 
as  follows: 

1.  In  no  cases  studied  was  the  shock  profile  found  to  be  steady 
in  time  as  is  generally  assumed  in  continuum  calculations.  ITie  cause 
of  the  nonsteady  behavior  is  the  propagation  of  well-defined  pulses 
(.solitary  waves)  behind  the  front  whose  proj)agation  velocity  varies 
with  amplitude,  i'he  difference  in  propagation  velocity  of  different 
solitons  introduces  a spreading  effect  which  apparently  prevents  the 
profile  from  ever  approaching  a steady  state. 

2.  For  the  initially  quiescent  lattice,  no  equilibration  of 
energy  occurs  behind  the  shock  front  and  there  is,  consequently,  no 
temperature  rise.  All  energy  introduced  by  Du*  shock  wave  is  either 
potential  or  ordered  translational  energy. 

3.  For  a lattice  with  a nonzero  initial  t emjierat ure , thermal i zat i on 
of  energy  behind  the  front  does  appear  to  occur  but  this  asiiect  of  the 
problem  requires  further  investigation.  Hie  profile,  however,  is  still 
nonsteady  and  the  transition  region  between  the  two  equilibrated  parts 

of  the  lattice  grows  with  time.  Again,  the  growth  results  from  the 
spreading  effect  of  solitary  waves  of  different  amplitudes  and  jiropa- 
gation  ve  loci  tie--. 

In  the  future,  we  intend  to  extend  the  calculations  to  a three- 
dimensional  model  to  see  if  similar  effects  will  persist  in  that  case. 

•1 0 

Recently,  Zabusky  ' h.is  interpreted  an  increase  in  thennal  conduction 
with  the  addition  of  anharmon i c i t les , observed  by  Payton,  Rich,  and 

V i sscher‘^^\  as  resulting  from  the  propagation  of  solitons  in  the  model 
studied.  That  model  was  two  dimensional  (as  well  as  isotopically 
impure)  and  it  therefore  appears  reasonable  to  expect  that  soli  ton 
propagation  might  be  an  imiiortant  effect  in  a three-dimensional  model 
as  well. 


4‘).  N..I.  Zalnisky,  "Solitons  and  Faierg)'  Iransport  in  Nonline;ir  Lattices", 

Computer  Phys.  i!omm.  S,  1 (i!)T.3). 

SO.  D.N.  Payton  111,  M.  Rich,  and  W.M.  Visscher,  "Lattice  niennal 
(iondiict  iv  i ty  in  Oisordered  ilarmonic  and  \nliarmonic  Crystal 
Models",  I’hys.  Rev.  KiO,  TOii  (l'i(i’). 


If  so,  we  believe  that  the  effects  may  be  significant  in  the  study 
of  shock-induced  detonations.  In  most  current  theories,  the  effect  of 
the  shock  front  is  ignored  entirely  and  the  only  function  the  shock  has 
is  to  raise  the  temperature,  pressure,  and  density  of  the  solid  to  higher- 
than-ambient  values.  Chemical  reactions  are  then  assumed  to  occur  in 
the  thermally  equilibrated  part  of  the  crystal  behind  the  front.  It 
would  appear,  however,  that  if  the  transition  region  grows  in  time, 
reactions  may  occur  in  a region  of  the  crystal  characterized  by  an 
extreme  state  of  nonequilibrium,  and  it  is  important  to  assess  the 
effects  of  this  nonequilibrium  environment  upon  chemical-reaction  rates. 

It  is  unlikely,  for  instance,  that  such  rates  can  be  represented  by  an 
Arrhenius-type  relation  whose  validity  is  dependent  upon  the  assumption 
of  equilibrium.  Furthermore,  the  solitons  at  the  shock  front  lead  to 
particle  velocities  significantly  higher  than  the  velocities  typically 
found  in  the  equilibrated  regions  of  the  compressed  crystal,  an  effect 
which  may  be  significant  in  understanding  the  initiation  process  in 
explosives . 

It  would  also  be  desirable  in  the  future  to  study  the  properties 
of  solitons  in  more  detail  than  has  been  possible  in  this  report.  Of 
particular  interest  would  be  to  study  the  effects  of  random  thermal 
oscillations  upon  the  propagation  of  solitons  of  different  amplitude. 

Such  a study  would  perhajis  lend  insight  into  why  (or  why  not)  equilibra- 
tion of  energy  occurs  behind  the  front  in  the  case  of  the  initially 
equilibrated  lattice.  It  would  also  be  interesting  to  study  the  effects 
of  perturbations  transverse  to  the  propagation  direction  of  the  soliton 
(in  a two-  or  three-dimensional  model).  It  would  appear  that  only  the 
absence  of  stability  under  such  perturbations  could  prevent  soliton 
propagation  from  being  an  important  effect  in  more  realistic,  three- 
dimensional  crystals. 
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